This article is accompanied by an Editorial View. Please see: Todd MM: EEGS, EEG processing, and the bispectral index. Anesthesiology 1998; 89:815-7.
Introduction-The Rationale for Monitoring
THE electroencephalogram (EEG) is enjoying a renaissance of interest as a clinical monitoring tool during anesthesia and sedation. This revival is the result of two recent events: first, retargeting the use of EEG from confirming deep surgical anesthesia to the assessment of lighter or sedative levels, and second, new technologic developments that have produced tangible progress in the creation of a monitor of “anesthetic depth.” This article is intended to review the technical bases of these developments as well as provide a synopsis of the relevant physiology involved in the generation and modulation of the EEG. A comprehensive examination of the clinical applications of EEG monitoring is, however, outside the scope of this review. A survey of EEG monitoring applications in anesthesiology has recently been published. [1]
The EEG was first described in 1875 by Richard Caton, [2]a physician in Liverpool, who noted electrical oscillations on the exposed cortical surface of animals. In 1929, Hans Berger, a psychiatrist in Jena, began a series of reports [3]that are commonly accepted as the first systematic description of human EEG. [4]Within 10 years Gibbs and Gibbs noted that the EEG was sensitive to presence of anesthetic agents. [5]The next 50 years brought significant improvement in the equipment that transduces, amplifies, and displays EEG. In the past 20 years some progress has been made in understanding the electrophysiology of the brain as it relates to the genesis of the EEG waveforms.
At present, EEG monitoring is used to assess, in real-time, either the state of “well-being” of the higher central nervous system (CNS) or the pharmacodynamic effect of an anesthetic drug. The EEG is widely accepted as being a highly sensitive, moderately specific indicator of CNS ischemia or hypoxia; thus EEG monitoring has been commonly used for this purpose during carotid surgery. [6-8]EEG monitoring for drug effect has three applications: a quantitative tool for the pharmacologic study of CNS-active agents;[9-11]the assessment of metabolic suppressive effect (e.g., dose control of thiopental for EEG burst suppression [12,13]); and recently, the assessment of CNS functional suppression (depth of sedation or anesthesia [14,15]).
The choice of physiologic end-point against which to correlate or test the EEG has proven critical. The first reported relationship compared EEG with changes in blood pressure after noxious stimulation. The degree of EEG depression before laryngoscopy correlated with magnitude of hemodynamic change after intubation. [16,17]However, later attempts to correlate EEG with gross purposeful movement in response to surgical incision have not uniformly reported positive results. [18-20]These inconsistent results may be the result of the anatomic and possible pharmacodynamic separation of the neural circuitry responsible for movement responses (spinal cord) from those responsible for the generation of the EEG signal (cerebrum). [21]EEG, as will be discussed, is a phenomenon of the rostral structures, particularly the cerebral cortex. Anesthetic-induced suppression of spinal function, i.e., surgical immobility, may be observable by monitoring spinal reflexes like F-waves. [22,23]When EEG is correlated with neural functions linked to the cortex such as awareness or memory, more reliable, clinically relevant results are obtained. [24,25]
The Genesis of the EEG
All bioelectric potentials observed on the skin are caused by the flow of ion-based electrical currents within the volume of the body. As in the case of the electrocardiogram (ECG), these macroscopic currents are the net summation of microscopic currents contributed by large populations of individual, electrically active cells.
The electrical activity of neurons may be divided into two categories: regenerative action potentials (AP) and postsynaptic potentials (PSP). PSPs occur when neurotransmitters released by a presynaptic neuron alter the permeability of ion channels in the postsynaptic neuron's cell membrane, altering its transmembrane ionic concentration gradients and thus its transmembrane voltage. The magnitude of an isolated PSP is proportional to the number of postsynaptic receptors that have bound agonist. Because neurotransmitter release is a localized phenomenon, the resulting changes in resting membrane potential also tend to be focal, with the magnitude of the voltage change decreasing exponentially with distance from the synapse and a membrane constant known as [Greek small letter lambda]. The “length constant”, [Greek small letter lambda], is analogous to a time constant in describing the exponential decay of a perturbation. In this case, [Greek small letter lambda] depends on the characteristics of the cell membrane and describes the distance along the membrane at which the voltage disturbance has decayed to 37% of Delta V0. The value of [Greek small letter lambda] is often in the range of 0.1-1.0 mm, thus rendering PSPs a localized phenomenon. The direction of the change in membrane potential can be either positive (depolarizing) or negative (hyperpolarizing) depending on which ionic species has its membrane permeability altered by the neurotransmitter. Synaptic activity thus creates focal patches of altered membrane potential, and ionic current flow occurs between these disturbances. PSPs slowly decay over time, bringing the membrane potential back to its resting value. The mechanism of decay is a combination of cessation of ligand-triggered channel activity, either resulting from removal of the neurotransmitter on inactivation of the ion channels or the PSP-induced currents that redistribute ionic charge to counteract the PSP. Decay times of PSPs range in 10s of milliseconds to seconds.
If the membrane potential of a neuron is depolarized beyond its intrinsic threshold value, an AP is initiated. APs propagate rapidly along the membrane without diminution in amplitude, sustained by voltage-sensitive sodium and potassium channels and the transmembrane concentration gradients of these species. Typically at a given point on the membrane, the excursion in potential caused by an AP lasts for about 2 ms and may reach approximately 100 mV in amplitude.
Cytoarchitecture
It is the relatively slowly changing currents from PSPs that are the basis of the macroscopic EEG. In isolation, the local current loops from a single neuron would be difficult to detect at the distance of the scalp. Fortuitously, the anatomy of the cerebral cortex provides a means of generating relatively robust signals. Cortical neurons are classified by their morphology. [26]A prevalent type is the pyramidal cell that has a long straight apical dendrite extending up through the cortical layers from the cell body directly toward the pial surface of its gyrus. Neighboring pyramidal cells therefore have roughly parallel dendrites. These pyramidal dendrites receive thousands of synaptic connections from other neurons, and neighboring dendrites tend to receive inputs from many of the same presynaptic sources and at homologous locations on the dendrite. Physical separation of inhibitory and excitatory PSPs on an individual apical dendrite creates bridging current loops between the PSPs, which are far larger than might be predicted from the length constants of dendrites (Figure 1). When neighboring pyramidal cells have similar and synchronous areas of altered membrane potentials, their current loops combine additively in the extracellular fluid to create a much larger regional current flow, one that can be detected by the voltage it creates on the scalp.
Figure 1. Formation of current loops in the extracellular space between pyramidal cell dendrites. Dendrites receive many simultaneous synaptically transmitted signals and will initiate an action potential when the net balance of the afferent excitatory traffic exceeds the inhibitory traffic by a margin that allows a localized patch of its membrane to depolarize to the critical value. If the excitatory activity is somewhat separated physically from the inhibitory activity, the resulting difference in membrane potential will cause a current to flow both inside and outside the membrane acting. The extracellular component of this current, when magnified by summation with the currents from neighboring parallel pyramidal dendrites, creates voltages that are detectable on the scalp as the EEG.
Figure 1. Formation of current loops in the extracellular space between pyramidal cell dendrites. Dendrites receive many simultaneous synaptically transmitted signals and will initiate an action potential when the net balance of the afferent excitatory traffic exceeds the inhibitory traffic by a margin that allows a localized patch of its membrane to depolarize to the critical value. If the excitatory activity is somewhat separated physically from the inhibitory activity, the resulting difference in membrane potential will cause a current to flow both inside and outside the membrane acting. The extracellular component of this current, when magnified by summation with the currents from neighboring parallel pyramidal dendrites, creates voltages that are detectable on the scalp as the EEG.
Control of Rhythm
As PSPs occur and decay, the EEG scalp voltage changes over time. Under ordinary circumstances, millions of PSPs are asynchronously firing all over the cortex, summing to create a complicated composite signal that cannot be decomposed back into component PSPs. Therefore, unlike the ECG, normal EEG has no obvious, repetitive patterns, nor does the shape of the EEG waveform correlate with specific underlying events. On the contrary, the EEG is a random-appearing signal. Decades of empirical observation indicate, however, that some statistical attributes derived from the EEG reflect and track the underlying state of the brain. Thus, there are some characteristics of the EEG that may be measured to provide a quantitative, if indirect, monitoring of some aspect of brain function. EEG monitoring for anesthesia-related purposes relies on this statistical approach. For example, the degree of PSP synchrony in pyramidal cells determines the net magnitude and frequency resulting from contributions of individual neurons. Higher cortical function is usually associated with desynchronization as neurons act more independently in the creation of conscious human behavior. Anesthesia and other mechanisms that depress consciousness are associated with increasing cortical synchrony. Anatomically, synchrony and level of consciousness are strongly influenced by neuronal circuit loops involving cortical connection with the brainstem and thalamus. [27]These circuits are sometimes called the EEG pacemaker.
Under some conditions the EEG may contain special stereotypic waveforms that can be used diagnostically. For example, “spikes” or “sharp waves” are sharply featured excursions in the EEG that are created by massive but usually transient synchrony. The presence of repetitive spikes is used in the diagnosis of epilepsy. Another example is burst suppression phenomenon (intermittent electrical activity interspersed with silence), which indicates a nonspecific (e.g., trauma, drugs, hypothermia, and so on) reduction in cerebral metabolic activity.
Once generated, the aggregate postsynaptic current flows must traverse the cerebrospinal fluid (CSF), the skull, and the scalp to be detected by skin surface electrodes. The CSF and scalp are relatively conductive compared with the skull, and the overall effect of transmission through these layers is a substantial spatial smearing of regional voltage differences. This means that the signal from an EEG electrode reflects activity over a wide area, not just directly under the electrode.
Signal Acquisition
Metal needle or gel electrodes are required to act as transducers, converting EEG (physiologic) ionic current into electronic current that may then be further processed by the EEG monitoring equipment.
Voltage signals, like the EEG, are always measured as a difference in potential between two points; thus a bioelectric amplifier has two signal inputs: a plus and a minus. Bioelectric amplifiers also have a third input for a reference electrode, which will be discussed. Because the electrical activity of the cortex is topographically heterogeneous, it is generally advantageous to measure this activity at several locations on the scalp. In diagnostic neurology, several systems of nomenclature for electrode placement have evolved; the most commonly used at present is the International 10-20 System. [28]The 10-20 system is based on meridians crossing the scalp based on key landmarks (the nasion, the inion, and the left and right aural tragus) with additional lines drawn over the midfrontal lobes and the midparietal lobes. The nomenclature uses a letter prefix designating brain site (i.e., C = central, F = frontal, P = parietal, T = temporal), and a number indicating the relative distance from the midline (nasion to inion), where right-sided electrodes are given even numbers and left-sided placements receive odd numbers. Electrodes on the midline are designated “Z”, i.e., FZ is a site overlying the falx between the frontal lobes, and P3 is a left parietal site.
Diagnostic EEG as performed in the Neurology clinic is seldom recorded with fewer than 16 channels (plus-minus pairs of electrodes) to localize abnormal activity. Monitoring 8 or 16 channels intraoperatively during carotid surgery is often recommended, although there is a paucity of data demonstrating increased sensitivity for the detection of cerebral ischemia compared with the two- or four-channel computerized systems more commonly available to anesthesia personnel. Although regional changes in EEG occur during anesthesia, [29]there is little evidence that these topographic features are useful markers of clinically important changes in anesthetic or sedation levels. [30]
Amplifiers
The EEG signal is but one of several voltage waveforms present on the scalp. In awake subjects, there are at least three other signals generated by physiologic processes: ECG (from the R-wave vector sweeping through the neck), electromyogram (EMG, from electromechanical activity of scalp muscle), and the electrooculogram (EOG, generated by movement of the electrical dipoles within the globes as the eyes move). In addition, the body acts as an antenna that picks up the powerline signal radiated by the cables in nearby walls and ceilings.
Although all these signals may contain interesting information,[dagger] if present, they distort and interfere with the EEG signal. An understanding of the essential characteristics of specific artifacts can be used to mitigate them. A well-designed bioelectric amplifier can remove or attenuate some of these signals as the first step in signal processing. The largest source of artifact is the powerline pickup. This artifact possesses two useful characteristics: it is the same over the entire body surface, and it is a single characteristic frequency. Because EEG voltage is measured as the potential difference between two electrodes placed on the scalp, both electrodes will have the same powerline artifact (i.e., it is a common mode signal). Common mode signals can be nearly eliminated in the electronics stage of an EEG machine by using a differential amplifier that has connections for three electrodes:“+,”“-,” and a reference. This type of amplifier detects two signals: the voltage between + and reference, and between - and reference, then subtracts the second signal from the first. The contribution of the reference electrode is common to both signals and thus cancels out. Attenuation of common mode artifact signals will be complete only if each the + and - electrodes are attached to the skin with identical contact impedance. If the electrodes do not have equal contact impedances, the amplitude of the common mode signal will differ, and they will not cancel exactly. Most commonly, the EEG is measured (indirectly) between two points on the scalp with a reference electrode on the ear or forehead. If the reference electrode is applied far from the scalp, i.e., the thorax or leg, there is always a chance that large common mode signals like the ECG will not be ideally canceled out, leaving some degree of a contaminating artifact.
Some artifacts, like the EMG, characteristically have most of their energy in a frequency range different from that of the EEG. Hence, the amplifier can band-pass filter the signal, passing the EEG and attenuating the EMG. Some EEG machines quantify and separately report EMG activity on the screen before filtering it from the EEG.
Signal Processing
Signal processing of an EEG is done to enhance and aid the recognition of some aspect of the EEG that correlates with the physiology and pharmacology of interest. Metaphorically, the goal is to separate this “needle” from an electrical haystack. The problem in EEG-based assessment of anesthetic state is that the characteristics of this needle are unknown, and because our fundamental knowledge of the CNS remains relatively limited, our needle-like constructs will, for the foreseeable future, be based on empirical observation. Assuming a useful quantitative EEG parameter (QEEG) is identified, it must be measured. The motivation for quantitation is threefold: to reduce the clinician's workload in analyzing intraoperative EEG, to reduce the level of specialized training to take advantage of EEG, and finally to develop a parameter that might, in the future, be used in an automated closed-loop titration of anesthetic or sedative drugs. This section will introduce some of the mechanics and mathematics behind signal processing.
Although it is possible to perform various types of signal enhancement on analog signals, the speed, flexibility, and economy of digital circuits have produced revolutionary changes in the field of signal processing. To use digital circuits, it is, however, necessary to translate an analog signal into its digital counterpart.
Analog signals are continuous and smooth. They can be measured or displayed with any degree of precision at any moment in time. The EEG is an analog signal; the scalp voltage varies smoothly over time.
Digital signals are fundamentally different because they represent discrete points in time and their values are quantitated to a fixed resolution rather than continuous. The binary world of computers and digital signal processors operates on binary numbers, which are sets of bits. A bit is quantal; it contains the smallest possible chunk of information: a single ON or OFF signal. More useful binary numbers are created by aggregating between 8 and 80 bits. The accuracy or resolution (q) of binary numbers is determined by the number of bits they contain: an 8-bit binary number can represent 28or 1 of 256 possible states at any given time; a 16-bit number 216or 65,536 possible states. If one were using an 8-bit number to represent an analog signal, the binary number would have, at best, a resolution of approximately 0.4%(1/256) over its range of measurement. Assuming, for example, the converter was designed to measure voltages in the range from -1.0 to +1.0 V, the step size of an 8-bit converter would be about 7.8 mV, and a 16-bit converter would be about about 30 [micro sign]V. EEG monitoring systems usually use 12-16 bits of resolution.
Digital signals are also quantized in time, unlike analog signals, which vary smoothly from moment to moment. When translation from analog to digital occurs, it occurs at specific points in time, and strictly speaking, the value of the resultant digital signal at all other points in time is indetermine (Figure 2). Translation from the analog to digital world is known as sampling or digitizing and in most applications it is set to occur at regular intervals. The reciprocal of the sampling interval is known as the sampling rate (fs) and is expressed in hertz (Hz or samples per second). A signal that has been digitized is commonly written as a function of a sample number, i, instead of analog time, t. For example, an analog voltage signal might be written at v(t), but after digitizing it would be referred to as v(i). Taken together, the set of sequential samples representing a finite block of time is referred to as an epoch or a realization. In statistical theory, the collection of all possible epochs produced from a given EEG state would be known as an ensemble.
Figure 2. Analog to digital conversion translates a smoothly continuous waveform into a Table offinite resolution numbers, x(k), representing the value of the original waveform at discrete points in time (in this example, the index k happens to equal integer seconds). This example illustrates sampling at 1.0 Hz with a resolution of four decimal digits.
Figure 2. Analog to digital conversion translates a smoothly continuous waveform into a Table offinite resolution numbers, x(k), representing the value of the original waveform at discrete points in time (in this example, the index k happens to equal integer seconds). This example illustrates sampling at 1.0 Hz with a resolution of four decimal digits.
The process of analog to digital translation inevitably leads to a loss of fidelity in the resulting digital signal. The resulting digital signal, x(i), can be thought of as the sum of (an impossibly) perfect digital copy of the true signal xu(i) plus an error term, e(i). The quantization error, e(i), is the difference between the sampled voltage and the true analog voltage. Quantization error can be reduced by increasing the number of bits used to represent the digitized sample. The signal processing designer must trade off increased accuracy against the increased cost of higher resolution hardware (including the A-to-D converter itself as well as a wider data path in the computing circuits (i.e., the computer arithmetic unit must be expanded to handle numbers with more bits) and more memory to hold the added bits.
When sampling is performed too infrequently, the fastest sine waves in the epoch will not be identified correctly. When this situation occurs, aliasing distorts the resulting digital data. Aliasing results from the requirement for a minimum of two points within a single cycle to identify a sinusoid. If sampling is not fast enough to place at least two sample points within a single cycle, the sampled wave will appear to be slower (longer cycle time) than the original. Aliasing is familiar to observers of the visual sampled-data system known as cinema. In a movie, where frames of a scene are captured at rate of approximate 24 Hz, rapidly moving objects like wagon wheel spokes often appear to rotate slowly or even backward. Aliasing in a digitized waveform is illustrated in Figure 3.
Figure 3. Aliasing is a distortion caused by an insufficient sampling rate that erroneously converts a fast analog wave into a slow digitized wave. This Figure beginswith a complex wave with a dominant frequency of approximately 0.43 Hz (i.e., approximately eight cycles in 20 s) and some smaller waves of even higher frequency. In this example, the signal is sampled at 0.25 Hz, which is less than one third the minimum rate needed to follow the changing waveform. The resulting digitized waveform, drawn with linear interpolation in between the sample points, has a frequency of approximately 0.05 Hz, clearly a poor representation of the original signal. Sampling at 1.0 Hz or faster would produce a digital waveform with useful fidelity. Note that linear interpolation (drawing a straight line between samples) used in the Figure isusually an unjustified assumption about how the analog signal changes between digitized samples.
Figure 3. Aliasing is a distortion caused by an insufficient sampling rate that erroneously converts a fast analog wave into a slow digitized wave. This Figure beginswith a complex wave with a dominant frequency of approximately 0.43 Hz (i.e., approximately eight cycles in 20 s) and some smaller waves of even higher frequency. In this example, the signal is sampled at 0.25 Hz, which is less than one third the minimum rate needed to follow the changing waveform. The resulting digitized waveform, drawn with linear interpolation in between the sample points, has a frequency of approximately 0.05 Hz, clearly a poor representation of the original signal. Sampling at 1.0 Hz or faster would produce a digital waveform with useful fidelity. Note that linear interpolation (drawing a straight line between samples) used in the Figure isusually an unjustified assumption about how the analog signal changes between digitized samples.
It is essential to always sample at a rate more than twice the highest expected frequency in the incoming signal (Shannon's sampling theorem [31]). Conservative design actually calls for sampling at a rate 4 - 10 times higher than the highest expected signal and to also use a low-pass filter before sampling to eliminate signals that have frequency components that are higher than expected. Low-pass filtering reduces high-frequency content in a signal, just like turning down the treble control on a stereo system. In monitoring work, EEG signals have long been considered to have a maximal frequency of 30 or 40 Hz, although 70 Hz is a more realistic limit. In addition, other signals present on the scalp include powerline interference at 60 Hz and the electromyogram, which, if present, may extend above 100 Hz. To prevent aliasing distortion in the EEG from these other signals, many digital EEG systems sample at a rate above 250 Hz (i.e., a digital sample every 4 ms).
Artifact Recognition
The problem of artifactual contamination must always be considered in EEG analysis. Artifact is particularly insidious in EEG analysis because even to the trained observer, much of true EEG resembles noise. Common artifacts include signals that have exceeded the dynamic range of the amplifier (voltage too big because of improper amplifier settings or movement of the electrodes on the skin). These artifacts are easy to recognize, but the epoch containing this artifact must be excluded from analysis because the original data cannot be reconstructed. Another common type of artifact, as noted previously, is caused by the presence of an additional signal that is outside the frequency range of the EEG. This type of signal might include electromyogram activity or powerline pickup. If the sampling rate is fast enough to avoid aliasing, these kinds of artifacts may be filtered out, leaving a still usable EEG signal. Other types of artifact, including the electrocardiogram and roller pump artifact (during cardiopulmonary bypass), occur within the frequency range of interest for EEG and may be recognized by their regularity. Anesthesia equipment, such as a train-of-four twitch stimulator or an evoked potential stimulator, may also create a patterned artifact in the EEG. In awake or lightly sedated subjects, eye blinks and rotation of the orbital globes create large, transient slow wave activity, which may be recognized based on the pattern of signal amplitude changes. A compendium of techniques for EEG artifact detection and mitigation is provided by Barlow. [32]In commercially available EEG monitors, epochs with artifact may be tagged but still processed or identified and excluded from further processing.
Time Domain Methods
The EEG is an alternating voltage composed of many wavelets (simple sine waves) superimposed on each other. Analysis of the EEG can be accomplished by examining how its voltage changes over time. This approach, known as time domain analysis, may use either a strict statistical calculation (i.e., the mean and variance of the sampled waveform or the median power frequency) or may use some ad hoc measurement based on the morphology of the waveform. Most of the commonly used time-domain methods are grounded in probabilistic analysis of “random” signals, and therefore some background on statistical approaches to signals is useful. Of necessity, the definitions of probability functions, expected values, and correlation are given mathematically and descriptively. However, the reader needs not feel compelled to attain a deep understanding of the equations presented here to continue. A more detailed review of the statistical approach to signal processing may be obtained from one of the standard texts. [33-35]Only one class of ad hoc time-domain methods, burst suppression quantitation, is in current use in perioperative monitoring systems, and it will be described below.
A few definitions related to the statistical approach to time-related data are required. The EEG is not a deterministic signal, which means that it is not possible to exactly predict future values of the EEG. Even if the exact future values of a signal cannot be predicted, some statistical characteristics of certain types of signals are predictable in a general sense. These roughly predictable signals are termed stochastic. The EEG is such a nondeterministic, stochastic signal because its future values can only be predicted in terms of a probability distribution of amplitudes already observed in the signal. This probability distribution, p(x), can be determined experimentally for a particular signal, x(t), by simply forming a histogram of all the observed values over time. A signal such as may be obtained by rolling dice has a probability distribution that is rectangular or uniform (i.e., the likelihood of all face values of a throw are equal, and in the case of a single die, p(x)= 1/6 for each possible value); a signal with a bell-shaped or normal probability distribution is termed gaussian. As illustrated in Figure 4, EEG amplitude histograms may have a nearly gaussian distribution. The concept of statistics, like the mean, standard deviation, skewness, and so on, used to describe a probability distribution will be familiar to many readers.
Figure 4. Probability density functions are simply histograms of the amplitude values versus the number of samples at each value in the sampled signal. The waveforms illustrated previously are four sequential 4-s epochs comprising 16 s of EEG in an awake subject. The voltage scale on the left is after amplification by a factor of 500,000. The signal was sampled at 128 Hz after anti-alias filtering with a bandpass of 0.3-30 Hz. The histograms on the right demonstrate that these examples of EEG are approximately gaussian, and because the probability curves are mildly dissimilar, the EEG is not strictly stationary. Nonetheless, this Figure illustratesthe point that different appearing waveform data may produce statistical descriptions that are relatively constant between epochs.
Figure 4. Probability density functions are simply histograms of the amplitude values versus the number of samples at each value in the sampled signal. The waveforms illustrated previously are four sequential 4-s epochs comprising 16 s of EEG in an awake subject. The voltage scale on the left is after amplification by a factor of 500,000. The signal was sampled at 128 Hz after anti-alias filtering with a bandpass of 0.3-30 Hz. The histograms on the right demonstrate that these examples of EEG are approximately gaussian, and because the probability curves are mildly dissimilar, the EEG is not strictly stationary. Nonetheless, this Figure illustratesthe point that different appearing waveform data may produce statistical descriptions that are relatively constant between epochs.
In statistical parlance, these measurements of a distribution belong to family of statistical parameters known as moments or expected values Ek. For any given probability distribution, p(x), a number of different descriptive moments can be calculated where each moment is identified by its order number, k. For an EEG signal (consisting of a sequential collection of digitized amplitude values x(i)), the family of moments is defined as:Equation 1
The first order moment (e.g., substituting k = 1 in Equation 1) is mathematically the arithmetic mean of signal ([micro sign]). Consider the example of the rolling die again; the mean, [micro sign], is simply the sum of each possible face value, x(i), each multiplied by its probability, p(x)= 1/6 in a fair die. Above the first order moment, it is preferable to calculate the central moment (moment of the signal after the mean value of the epoch, [micro sign], has been subtracted out), as defined in Equation 10.
The second order central moment (k = 2) is the variance (i.e., the square of the standard deviation), the third central moment (k = 3) is the skewness, and the fourth is the kurtosis of the probability function.
If the probability function, p(x), of a stochastic signal, x(i), does not change over time, that process is stationary. The EEG is not strictly stationary as its statistical moments may change significantly within seconds (Figure 4), or it may be stable for 10s of minutes (quasistationary). [36,37]If the EEG is at least quasistationary, then it may be reasonable to check it for the presence of rhythmicity, where rhythmicity is defined as repetition of patterns in the signal. Patterns can be identified mathematically using the concept of correlation. Correlation between two signals measures the likelihood of change in one signal leading to a consistent change in the other. In assessing the presence of rhythms, autocorrelation is used, which tests the match of the original signal against different starting time points of the same signal. If rhythm is present, then at a particular offset time (equal to the interval of the rhythm), the correlation statistic increases, suggesting a repetition of the original signal voltage. The autocorrelation of signal x (i.e., correlation of x vs. x) is denoted as [Greek small letter gamma]xx([Greek small letter tau]) where [Greek small letter tau] is the offset time interval or lag.
Empirically it is known that the EEG has a mean voltage of zero, over time: it is positive as often as it is negative. However, the EEG and its derived statistical measurements seldom have a true gaussian probability distribution. This observation complicates the task of a researcher or of some future automated EEG alarm system that seeks to identify changes in EEG over time. Strictly speaking, non-gaussian signals should not be compared using the common statistical tests, such as t tests or analysis of variance, that are appropriate for normally distributed data. Instead, there are three options: non-parametric statistical tests, a transform to convert non-gaussian EEG data to a normal (guassian) distribution, or higher order spectral statistics (see below). Transforming non-gaussian data by taking their logarithm is frequently all that is required to allow analysis of the EEG as a normal distribution. [38]For example, a brain ischemia detection system may try to identify when slow-wave activity has significantly increased. A variable like “delta” power (described below), which measures slow-wave activity, has a highly non-gaussian distribution; thus, directly comparing this activity at different times requires the non-parametric Kruskal-Wallis or Friedman's test. However, the logarithm of delta power may produce a nearly normal p(x) curve; therefore the more powerful analysis of variance with repeated measures could be used appropriately to detect changes in log(delta power) over time. Log transformation is not a panacea, however, and whenever statistical comparisons of QEEG are to be made, the data should be examined to verify the assumption of normal distribution.
Clinical Applications of Time Domain Methods
Historically the first intraoperative application of EEG analysis used time domain-based methods. In 1950, Falconer and Bickford noted that the electrical power in the EEG (power = voltage [middle dot] current = voltage (2/resistance)) was associated with changes in the rate of thiopental or ethyl ether administration. Using analog technology, they computed a power parameter as (essentially) a moving average of the square of EEG voltage and used it to control the flow of diethyl ether to a vaporizer. This system was reported to successfully control depth of anesthesia in 50 patients undergoing laparotomy. [39]Digital total power (TP = sum of the squared values of all the EEG samples in an epoch) was later used by several investigators, but it is known have several problems, including its sensitivity to electrode location and its insensitivity to important changes in frequency distribution. Arom [40]reported that a decrease in TP may predict neurologic injury after cardiac surgery.
Hjorth [41]created a trio of combinations of conventional (time domain-based) descriptive statistics parameters: activity, mobility, and complexity. Activity is defined as the variance of the signal amplitude within an epoch, i.e., the square of the standard deviation of the digitized data points, and provides a measure of the mean power of the signal. Mobility can be considered an approximation of the average frequency of the EEG. It is defined as the standard deviation of the first derivative of the EEG signal (i.e., the intersample slope of the waveform) divided by the standard deviation of the original signal. Complexity is a variable to quantify the degree of curve complexity beyond a baseline sine wave. These parameters of Hjorth are frequently used in the related application of sleep staging [42,43]but have not directly been tested in perioperative monitoring.
A time domain-based approach to analysis of the frequency information within the EEG was reported by Burch [44]and by Klein and Davis, [45]who estimated an “average” frequency by detecting the number of times the EEG voltage crosses the zero voltage level per second. Investigators have not reported strong clinical correlations with zero crossing frequency (ZXF). Although simple to calculate in the era before inexpensive computer chips, the ZXF parameter is not simply related to frequency domain estimates of frequency content as demonstrated in Figure 5because not all waves in the signal will cross the zero point. Demetrescu [46]refined the zero crossing concept to produce what he termed aperiodic analysis. This method simply splits the EEG into two frequency bands (0.5 - 7.9, 8 - 29.9 Hz), and the filtered waveforms from the high and low frequency bands are each separately sent to a relative minima detector. Here, a wavelet is defined as a voltage fluctuation between adjacent minima, and its frequency is defined as the reciprocal of the time between the waves. Wavelet amplitude is defined as difference between the intervening maxima and the average of the two minima voltages. Aperiodic analysis produces a spectrum-like display called a glass box (Figure 6), which plots a sampling of detected wavelets as an array of “telephone poles” whose height represents measured wave amplitude, distance from the left edge frequency (in a logarithmic scale), and distance from the lower edge time since occurrence. The Lifescan Monitor (Diatek, San Diego, CA) was an implementation of aperiodic analysis; it is not commercially available at present, but the algorithms are described in detail in an article by Gregory and Pettus. [46]
Figure 5. The zero crossing algorithm measures the intervals between sequential crossings of the zero voltage axis. This technique is simple to implement and is an accurate measure of mean frequency when the EEG signal consists of waves that are similar in frequency and amplitude to each other, as in part A. The zero crossing algorithm has poor accuracy when the EEG contains a range of wave amplitudes and frequencies. As illustrated in part B, small fast waves riding on a larger slow wave do not cross the zero axis and are ignored.
Figure 5. The zero crossing algorithm measures the intervals between sequential crossings of the zero voltage axis. This technique is simple to implement and is an accurate measure of mean frequency when the EEG signal consists of waves that are similar in frequency and amplitude to each other, as in part A. The zero crossing algorithm has poor accuracy when the EEG contains a range of wave amplitudes and frequencies. As illustrated in part B, small fast waves riding on a larger slow wave do not cross the zero axis and are ignored.
Figure 6. Aperiodic display as implemented on the Lifescan [registered sign](Diatek, San Diego, CA). This time domain display depicts a 4-min trend where the newest information is drawn at the bottom of the near face of the “glass box.” Wavelets are depicted as vertical bars whose height is proportional to the wavelet amplitude and whose position along the x axis is proportional to its frequency band (roughly log10(frequency) horizontally) and recedes in depth with passing time. Frequency band definitions are given in Table 1. In this example, there was activity in both the delta and alpha bands until 2 min before the current time, when the alpha activity ceased.
Figure 6. Aperiodic display as implemented on the Lifescan [registered sign](Diatek, San Diego, CA). This time domain display depicts a 4-min trend where the newest information is drawn at the bottom of the near face of the “glass box.” Wavelets are depicted as vertical bars whose height is proportional to the wavelet amplitude and whose position along the x axis is proportional to its frequency band (roughly log10(frequency) horizontally) and recedes in depth with passing time. Frequency band definitions are given in Table 1. In this example, there was activity in both the delta and alpha bands until 2 min before the current time, when the alpha activity ceased.
Burst Suppression Quantitation
During deep anesthesia, the EEG may develop a peculiar pattern of activity, which is evident in the time domain signal. This pattern, known as burst suppression, is characterized by alternating periods of normal to high voltage activity changing to low voltage or even isoelectricity rendering the EEG inactive in appearance. After head trauma or brain ischemia, this pattern carries a grave prognosis; however, it may also be induced by large doses of general anesthetics, in which case, burst suppression has been associated with reduced cerebral metabolic demand and possible brain “protection” from ischemia. Titration to a specific degree of burst suppression has been recommended as a surrogate endpoint against which to titrate barbiturate coma therapy. The burst suppression ratio (BSR) is a time domain EEG parameter developed to quantitate this phenomenon. [47,48]To calculate this parameter, suppression is recognized as those periods longer than 0.50 s, during which the EEG voltage does not exceed approximately +/− 5.0 mV. The time in a suppressed state is measured, and the BSR is report as the fraction of the epoch length where the EEG is suppressed (Figure 7).
Figure 7. The burst suppression ratio algorithm is a time domain analysis technique that quantifies the degree of burst suppression. In this case, the threshold for suppressions is +/− 5.0 [micro sign]V. The BSR is calculated as the sum of intervals of suppression each at least 0.5 s long divided by the epoch length.
Figure 7. The burst suppression ratio algorithm is a time domain analysis technique that quantifies the degree of burst suppression. In this case, the threshold for suppressions is +/− 5.0 [micro sign]V. The BSR is calculated as the sum of intervals of suppression each at least 0.5 s long divided by the epoch length.
The random character of the EEG dictates that QEEG parameters extracted will exhibit a moment-to-moment variation without discernible change in the patient's state. Thus, output parameters are often smoothed by a moving average before display. Because of the particularly variable (non-stationary) nature of burst suppression, the BSR should be averaged over at least 15 epochs (60 s).
Frequency Domain Methods
An important alternative approach to time domain analysis examines signal activity as a function of frequency. So-called frequency domain analysis has evolved from the study of simple sine waves. A simple sine (or cosine) wave is a function of time, t, that can be completely described by Equation 2with three parameters: amplitude, frequency, and phase angle.
Amplitude is one half the peak-to-peak voltage; frequency is number of complete cycles per second; and phase angle is the way to describe the starting point of the waveform. The phase angle can be understood by considering a sine wave as the vertical displacement of a spoke on a spinning wheel, as illustrated in Figure 8. The convention is that the spoke starts in the horizontal position (0 [degree sign]), starting the sine at the midpoint of its vertical excursion. If upright, the spoke would be considered at 90 [degree sign] phase. In this model of a sine wave generator, the radius of the wheel is the sine amplitude, and the number of revolutions per second is the sine wave frequency.
Figure 8. A rotating vector or spoke describes a sinusoid over time. In this example, image a wheel, rotating counterclockwise, with a light source in its rim adjacent to the marked spoke. As the wheel turns, a graph of the vertical position of the light versus time will produce the indicated sine wave. The rotational speed of the wheel determines the frequency of the resulting sine wave. The initial angle of the spoke is the phase angle of the sine wave as used in Equation 2. In this illustration, the wheel starts at three different phase angles. Note that the sine frequency is independent of the phase angle.
Figure 8. A rotating vector or spoke describes a sinusoid over time. In this example, image a wheel, rotating counterclockwise, with a light source in its rim adjacent to the marked spoke. As the wheel turns, a graph of the vertical position of the light versus time will produce the indicated sine wave. The rotational speed of the wheel determines the frequency of the resulting sine wave. The initial angle of the spoke is the phase angle of the sine wave as used in Equation 2. In this illustration, the wheel starts at three different phase angles. Note that the sine frequency is independent of the phase angle.
Returning to EEG signals, it can be demonstrated (Fourier's theorem) that any arbitrarily complicated time-varying waveshape, x(t), can be decomposed into the sum of simple sine or cosine waves. Formally this is true as long as the original waveform was repetitive (e.g., x(t +[Greek small letter tau])= x(t), where [Greek small letter tau] is the period of repetition) and of infinite length. The value of [Greek small letter tau] is important because it determines the fundamental frequency, (f0= 2 [Greek small letter pi]/[Greek small letter tau]), of the signal. Informally, nonrepetitive signals are commonly decomposed with little error if the following caveats are noted. In practice, waveforms are finite and may not be truly repetitive. Nonperiodic signals can be treated as though the period of repetition has become infinitely long ([Greek small letter tau][arrow right][infinity]). The fundamental frequency then becomes infinitely small as does the separation of spectral components; thus the spectrum of nonperiodic signal is a continuous function.
The Fourier Equation (Equation 3) states that a periodic signal, x(t), consists of discrete components: a DC (zero frequency) amplitude, A0, plus sinusoids at the fundamental frequency, f0. and at an infinite number of integer multiples of the fundamental frequency, mf0(i.e., harmonics).
This rule applies whether the signal is in an analog or digital form. Figure 9demonstrates this process with the synthesis of a square wave from sinusoidal harmonics. The series of coefficients, Amand Bm, comprise the amplitude spectrum of the original signal, x(t), and the [Greek small letter theta]mcoefficients form the phase spectrum.
Figure 9. An example of the Fourier theorem that any repetitive waveform can be constructed from, or deconstructed to, a series of simple sine wave harmonics. Harmonics are sine waves whose frequencies are integer multiples of the slowest component wave. In this case, a square wave is constructed by sequentially adding odd harmonics. As more harmonics are added, the reconstructed wave becomes closer to an ideal square wave. The final square wave, x(t), can be written as:Equation 11where f0is the fundamental frequency of the square wave and a0is its amplitude. In the case of a square wave, all the harmonic phase angles are zero.
Figure 9. An example of the Fourier theorem that any repetitive waveform can be constructed from, or deconstructed to, a series of simple sine wave harmonics. Harmonics are sine waves whose frequencies are integer multiples of the slowest component wave. In this case, a square wave is constructed by sequentially adding odd harmonics. As more harmonics are added, the reconstructed wave becomes closer to an ideal square wave. The final square wave, x(t), can be written as:Equation 11where f0is the fundamental frequency of the square wave and a0is its amplitude. In the case of a square wave, all the harmonic phase angles are zero.
A Fourier analysis generates a frequency spectrum, which is simply a histogram of amplitudes or phase angles as a function of frequency. The concept is well illustrated by the effect of passing a white light through a glass prism, creating a rainbow (or spectrum). Each color of light represents a unique frequency photon, and the relative brightness among the colors is a measure of the energy amplitude at each frequency. Any measured signal transformed by the Fourier technique into the frequency domain will have both an amplitude and a phase component for each harmonic frequency. Mathematically it is usually expedient to describe both components with a single complex number.
To review, regular (real) numbers can be considered as describing a point on a number line. A new number line drawn at right angles to the real number line defines the “imaginary” number line, and the plane defined by the intersecting number lines is known as “complex” number space. A complex number, z, is a vector describing a point on the complex plane. The vector has two components, x, which is the projection of the vector on the real number axis, and jy, which is the projection on the imaginary axis (j is the right angle vector for an imaginary number, and it is defined as the square root of -1). An operation frequently performed in signal processing is complex conjugation where the sign of the imaginary component is inverted, i.e., if z = x + jy, then the conjugate z*= x - jy. The complex number vector can also be described as a magnitude (i.e., vector length [vertical bar] z [vertical bar]=[radical][chi squared]+ y2) and a phase (angle subtended to the real number axis [Greek small letter phi]= atan). This complex vector format provides the basis for interchangeability between complex numbers and the “spinning” wheel model as described previously for the phases of sinusoids. The Fourier equation (Equation 3) is thus written using the compact notation of complex exponentials:Equation 4
In Equation 4, m, again is the harmonic number, and cmis a complex Fourier coefficient that contains amplitude and phase information for each harmonic.
The conversion of a time domain waveform, x(t), into its sine was frequency components, X(f), is known as a Fourier transformation. This transformation, under ideal conditions, does not alter or reduce the information content within the waveform, and an inverse Fourier transformation will reconstitute the original waveform (i.e., the transformation is symmetric). Performing the Fourier transformation is then just computing the series of amplitude coefficients, Amand Bm, and the harmonic phase angles, [Greek small letter theta]m, from Equation 3, or the series of complex coefficients, cmper Equation 4.
The complex number version of the Fourier equation (Equation 4) defines a time domain signal to be a sum of complex harmonics. One may easily separate the amplitude components of the spectrum by taking the magnitude of each of the complex magnitudes, [vertical bar] cm[vertical bar], which provides a set of values that comprise the amplitude spectrum. Squaring the values of amplitude spectrum creates the power spectrum. The phase spectrum is the set of [Greek small letter phi]mcoefficients. Recall that m is the harmonic, a function of frequency expressed as multiple of the fundamental frequency.
The conversion of an analog signal to a sampled digital signal further modifies the signal spectrum. If a set of N digital samples from an analog waveform creates an epoch that is one complete period of a periodic waveform, then the discrete Fourier transform (DFT) spectrum will consists of N spectral components:Equation 5where n = 0, 1, 2,…, N-1 for all harmonics, and k is the sample index.
(Equation 5) is the digital version of the Fourier transform and uses the signal processing convention that time domain signals are written in lower case, x(t) or x(kT), and the equivalent complex frequency domain spectrum in upper case, X(f).
Increasing the duration of the sampling epoch increases the number of digital samples of time domain data, which in turn, increases the number of resulting frequency domain spectral components. The frequency resolution of the DFT is the reciprocal of the epoch duration, e.g., the Fourier spectrum produced by a 4-s epoch has a frequency resolution of 0.25 Hz.
The arbitrary slicing of a continuous signal stream into finite epochs (an epoch consists of a sequential set of samples) can introduce contamination in the form of artifactual frequencies created by the abrupt transitions at the ends of the epoch. For example, a digitized epoch of a sine wave starting and ending at a phase angle of 0 [degree sign] would have zero-valued samples at either end of the epoch, and the computed spectrum would be an accurate representation of the true frequency content of the signal, in this case, a single non-zero value at the frequency of the sine wave. In the case of a sampled sine wave starting and ending at 90 [degree sign] phase, the sampled values abruptly start and end at a non-zero value. The non-zero start and finish points cause the frequency domain to artifactually add high frequency components to the spectrum to match the abrupt, step-like transitions in the time domain signal.
This type of distortion is minimized by multiplying each time domain amplitude point within the epoch against the corresponding value of a window function. A window function is a numerical series containing the same number of elements as the number of signal samples in the epoch. Window element values tend toward zero at both ends and toward unity in the middle. A variety of window functions have been described, including the rectangle, the triangle, the Hanning, and the Blackman functions. [49]The effect of windowing is illustrated in Figure 10A. Until recently, commercial EEG analyzers did not perform windowing on sampled data. An illustration of Fourier transform-based EEG analysis is provided in Figure 10B.
Figure 10. (A) An illustration of time domain-based EEG processing. The top waveform is the original signal (after analog anti-alias filtering with a bandpass of 0.3-30 Hz) digitized at 256 Hz. The middle tracing demonstrates the effect of windowing on the original signal. Windowing is a technique used to reduce distortion from epoch end artifacts in subsequent frequency domain processing. A window consists of a set of digital values with the same number of members as the data epoch. In this case, a Blackman window was used. The Blackman window is defined as:Equation 12where i is the sample number and n is the number of samples in the epoch. The window operation multiplies each data sample value against its corresponding window value, i.e., the resulting waveform z(i)= x(i)[middle dot] w(i) for each value of i in the epoch. The bottom tracing is the autocorrelation function of this epoch of EEG. The autocorrelation provides much of the same information as a frequency spectrum because it can identify rhythmicities in the data. In this case, the strongest autocorrelation is at time = 0 as might be anticipated, and there are some weak rhythmicities that taper off as the lag increases above 1 s. (B) Continuing with the same epoch of digitized EEG, the top two tracings are the real and imaginary component spectra, respectively, resulting from the Fourier transform. The middle trace is the phase spectrum, which is usually discarded because of the present lack of known clinically useful correlation. The bottom tracing is the power spectrum. It is calculated as the sum of the squared real and imaginary components at each frequency (i.e., measuring the squared magnitude for each frequency value of the complex spectrum, X(f). Recall that power equals squared voltage. Note that the power spectrum, by reflecting only spectral magnitude, has explicitly removed whatever phase versus frequency information present in the original complex spectrum. From the power spectrum, the QEEG and relative band powers are calculated as described in the text and in Table 1.
Figure 10. (A) An illustration of time domain-based EEG processing. The top waveform is the original signal (after analog anti-alias filtering with a bandpass of 0.3-30 Hz) digitized at 256 Hz. The middle tracing demonstrates the effect of windowing on the original signal. Windowing is a technique used to reduce distortion from epoch end artifacts in subsequent frequency domain processing. A window consists of a set of digital values with the same number of members as the data epoch. In this case, a Blackman window was used. The Blackman window is defined as:Equation 12where i is the sample number and n is the number of samples in the epoch. The window operation multiplies each data sample value against its corresponding window value, i.e., the resulting waveform z(i)= x(i)[middle dot] w(i) for each value of i in the epoch. The bottom tracing is the autocorrelation function of this epoch of EEG. The autocorrelation provides much of the same information as a frequency spectrum because it can identify rhythmicities in the data. In this case, the strongest autocorrelation is at time = 0 as might be anticipated, and there are some weak rhythmicities that taper off as the lag increases above 1 s. (B) Continuing with the same epoch of digitized EEG, the top two tracings are the real and imaginary component spectra, respectively, resulting from the Fourier transform. The middle trace is the phase spectrum, which is usually discarded because of the present lack of known clinically useful correlation. The bottom tracing is the power spectrum. It is calculated as the sum of the squared real and imaginary components at each frequency (i.e., measuring the squared magnitude for each frequency value of the complex spectrum, X(f). Recall that power equals squared voltage. Note that the power spectrum, by reflecting only spectral magnitude, has explicitly removed whatever phase versus frequency information present in the original complex spectrum. From the power spectrum, the QEEG and relative band powers are calculated as described in the text and in Table 1.
Fast Fourier Transform
The original integral-based approach to computing a Fourier transform is computationally laborious, even for a computer. In 1965, Cooley and Tukey published an algorithm for efficient computation of Fourier series from digitized data. [50]This algorithm is known as the fast Fourier transform (FFT). More information about the implementation of FFT algorithms may be found in the text by Brigham. [51]Whereas the DFT of a sequence of N data points requires N2complex multiplications (a relatively time-consuming operation for a microprocessor), the FFT requires only N [middle dot](log2N)/2 complex multiplications. When the number of points is large, the difference in computation time is significant; for example, if N = 1,024, the FFT is faster by a factor of about 200.
In clinical monitoring applications, the results of a EEG Fourier transform are graphically displayed as a power versus frequency histogram, and the phase spectrum has been traditionally discarded as uninteresting. Whereas the frequency spectrum is relatively independent of the start point of an epoch (relative to the waveforms contained), the Fourier phase spectrum is highly dependent on the start point of sampling and thus very variable. Spectral array data from sequential epochs are plotted together in stack (like pancakes) so that changes in frequency distribution over time are readily apparent. Raw EEG waveforms, because they are stochastic, cannot be usefully stacked together because the results would be a random superposition of waves. However, the EEG's quasistationarity creates spectral data that are relatively consistent from epoch to epoch, allowing enormous visual compression of spectral data by stacking and thus simplified recognition of time-related changes in the EEG. Consider that raw EEG is usually plotted at a rate of 30 mm/s or 300 pages per hour, whereas the same hour of EEG plotted as a spectral array could be examined on a single page.
There are two types of spectral array displays available in commercial instruments: the compressed spectral array (CSA) and the density spectral array (DSA). The CSA presents the array of power versus frequency versus time data as a pseudo-three-dimensional topographic perspective plot (Figure 11), and the DSA presents the same data as a gray scale-shaded or colored two-dimensional contour plot. Although both convey the same information, the DSA is more compact, whereas the CSA permits better resolution of the power or amplitude data.
Figure 11. The creation of a spectral array display involves the transformation of time domain raw EEG signal into the frequency domain via the fast Fourier transform. The resulting spectral histograms are smoothed and plotted in perspective with hidden-line suppression for CSA displays (on left) or by converting each histogram value into a gray value for the creation of a DSA display (on right).
Figure 11. The creation of a spectral array display involves the transformation of time domain raw EEG signal into the frequency domain via the fast Fourier transform. The resulting spectral histograms are smoothed and plotted in perspective with hidden-line suppression for CSA displays (on left) or by converting each histogram value into a gray value for the creation of a DSA display (on right).
Bispectrum
The effort to glean useful information from the EEG has led from first order (mean and variance of the amplitude of the signal waveform) to second order (power spectrum, or its time domain analog, autocorrelation) statistics and now to higher order statistics. Higher order statistics include the bispectrum and trispectrum (third and fourth order statistics, respectively). Little work has been published to date on trispectral applications in biology, but there have been more than 100 papers and abstracts to date related to bispectral analysis of the EEG. Whereas the phase spectrum produced by Fourier analysis measures the phase of component frequencies relative to the start of the epoch, the bispectrum measures the correlation of phase between different frequency components as described below. What exactly these phase relationships mean physiologically is uncertain; one simplistic model holds that strong phase relationships relate inversely to the number of independent EEG pacemaker elements. Bispectral analysis has several additional characteristics that may be advantageous for processing EEG signals: gaussian sources of noise are suppressed, thus enhancing the signal-to-noise ratio for the non-gaussian EEG, and bispectral analysis can identify non-linearities, which may be important in the signal generation process. A complete treatment of higher order spectra may be found in the text by Proakis et al. [52]
As noted previously, the bispectrum quantifies the relationship among the underlying sinusoidal components of the EEG. Specifically, bispectral analysis examines the relationship between the sinusoids at two primary frequencies, f1and f2, and a modulation component at the frequency f1+ f2. This set of three frequency components is known as a triplet (f1, f2, and f1+ f2). For each triplet, the bispectrum, B(f1, f2), a quantity incorporating both phase and power information, can be calculated as described below. The bispectrum can be decomposed to separate out the phase information as the bicoherence, BIC(f (1), f2), and the joint magnitude of the members of the triplet, as the real triple product, RTP(f1, f2). The defining equations for bispectral analysis are described in detail below.
A high bicoherence value at (f1, f2) indicates that there is a phase coupling within the triplet of frequencies f1, f2, and f1+ f2. Strong phase coupling implies that the sinusoidal components at f1and f2may have a common generator, or that the neural circuitry they drive may, although some nonlinear interaction, synthesize a new, dependent component at the modulation frequency, f1+ f2. An example of such phase relationships and the bispectrum is illustrated in Figure 12.
Figure 12. The bispectrum is calculated in a two- dimensional space of frequency1versus frequency2as represented by the coarsely cross-hatched area. Because of the symmetric redundancy noted in the text and the limit imposed by the sampling rate, the bispectrum need only be calculated for the limited subset of frequency combinations illustrated by the darkly shaded triangular wedge. A strong phase relationship between f1, f2, and f1+2creates a large bispectral value B(f1, f2) represented as a vertical spike rising out of the frequency versus frequency plane. In panel A, three waves having no phase relationship are mixed together producing the waveform shown in the upper right. The bispectrum of this signal is everywhere equal to zero. In panel B, two independent waves at 3 and 10 Hz are combined in a non-linear fashion, creating a new waveform that contains the sum of the originals plus a wave at 13 Hz, which is phase-locked to the 3- and 10-Hz components. In this case, computation of the bispectrum reveals a point of high bispectral energy at f (1)= 3 and f2= 10 Hz.
Figure 12. The bispectrum is calculated in a two- dimensional space of frequency1versus frequency2as represented by the coarsely cross-hatched area. Because of the symmetric redundancy noted in the text and the limit imposed by the sampling rate, the bispectrum need only be calculated for the limited subset of frequency combinations illustrated by the darkly shaded triangular wedge. A strong phase relationship between f1, f2, and f1+2creates a large bispectral value B(f1, f2) represented as a vertical spike rising out of the frequency versus frequency plane. In panel A, three waves having no phase relationship are mixed together producing the waveform shown in the upper right. The bispectrum of this signal is everywhere equal to zero. In panel B, two independent waves at 3 and 10 Hz are combined in a non-linear fashion, creating a new waveform that contains the sum of the originals plus a wave at 13 Hz, which is phase-locked to the 3- and 10-Hz components. In this case, computation of the bispectrum reveals a point of high bispectral energy at f (1)= 3 and f2= 10 Hz.
Calculation of the bispectrum of a digitized epoch, x(i), begins with an FFT to generate complex spectral values, X(f). For each possible triplet, the complex conjugate of the spectral value at the modulation frequency X*(f1+ f2) is multiplied against the spectral values of the primary frequencies of the triplet (Equation 6). This multiplication is the heart of the bispectral determination: if at each frequency in the triplet, there is a large spectral amplitude (i.e., a sinusoid exists for that frequency) and if the phase angles for each are aligned, then the resulting product will be large; if one of the component sinusoids is small or absent or if the phase angles are not aligned, the product will be small. Finally, the complex bispectrum is converted to a real number by computing the magnitude of the complex product. Equation 6
If one starts by sampling EEG at 128 Hz into 4-s epochs, then the resulting Fourier spectrum will extend from 0 to 64 Hz at 0.25 Hz resolution, or a total of 256 spectral points. If all triplets were to be calculated, there would be 65,536 (256 [middle dot] 256) points. Fortunately, it is unnecessary to calculate the bispectrum for all possible frequency combinations. The unique set of frequency combinations to calculate a bispectrum can be visualized as a wedge of frequency versus frequency space (Figure 12). The combinations outside this wedge need not be calculated because of symmetry (i.e., B(f1, f2)= B(f2, f1)) and because the range of allowable modulation frequencies, f1+ f2, is limited to frequencies less than or equal to half the sampling rate. Still, because this calculation must be performed, using complex number arithmetic, for at least several thousand triplets, it is easy to see that it is a major computational burden.
As noted previously, computation of the bispectrum itself is only the beginning for complete higher order spectral analysis. If one is interested in isolating and examining solely the phase relationships, as noted previously, the bispectrum must have the existing variations in signal amplitude normalized. Recall that the amplitude of a particular Fourier spectral element X(f) is determined by the magnitude or the length of its complex number vector. The RTP (Equation 7) is formed from the multiplied product of the squared magnitudes of the three spectral values in the triplet:
The square root of the RTP yields the joint amplitude of the triplet, the factor that is used to normalize the bispectrum into the bicoherence. The bicoherence, BIC(f1, f2), (Equation 8) is a number that varies from 0 - 1 in proportion to the degree of phase coupling in the frequency triplet:
(Figure 13) illustrates some representative data during bispectral analysis. Computing the bispectrum of a stochastic biological signal, like the EEG, generally requires that the signal be divided into relatively short epochs for calculation of bispectrum and bicoherence, which are then averaged over a number of epochs to provide a relatively stable estimate of the true bispectral values.
Figure 13. Examples of bispectral analysis. (A) An epoch of raw EEG that has been anti-alias filtered and sampled at 125 Hz. (B) The power spectrum resulting from Fourier analysis of the epoch in part A. In this epoch, the large contribution from frequencies around 1 Hz overwhelms the contributions from higher frequencies. (C) The same data as part B are plotted with a logarithmic power scale, allowing the contribution from faster waves to be identified despite the presence of large slow waves. (D) After computing the complex spectrum X(f) for this epoch of data, the bispectrum may be calculated. As described in Equation 5, the bispectrum is computed as a function of two frequencies f1and f2. Part D of this Figure illustratesthe resulting two-dimensional contour plot of the bispectrum of the same epoch. The largest bispectral features here are in the upper left corner, relating pairs of delta-range frequency pairs. The amplitude scale is again logarithmic with the heavy lines indicating powers of 10. (E) The same bispectral data as part C plotted in three-dimensional perspective. The frequency scales are the same as part C. A plot of the real triple product function does not appear to be significantly different than this plot. The amplitude scale here and in panel E is linear. (F) However, when the bicoherence, BIC, is plotted here, phase coherence features are visible in a wide distribution of frequency.
Figure 13. Examples of bispectral analysis. (A) An epoch of raw EEG that has been anti-alias filtered and sampled at 125 Hz. (B) The power spectrum resulting from Fourier analysis of the epoch in part A. In this epoch, the large contribution from frequencies around 1 Hz overwhelms the contributions from higher frequencies. (C) The same data as part B are plotted with a logarithmic power scale, allowing the contribution from faster waves to be identified despite the presence of large slow waves. (D) After computing the complex spectrum X(f) for this epoch of data, the bispectrum may be calculated. As described in Equation 5, the bispectrum is computed as a function of two frequencies f1and f2. Part D of this Figure illustratesthe resulting two-dimensional contour plot of the bispectrum of the same epoch. The largest bispectral features here are in the upper left corner, relating pairs of delta-range frequency pairs. The amplitude scale is again logarithmic with the heavy lines indicating powers of 10. (E) The same bispectral data as part C plotted in three-dimensional perspective. The frequency scales are the same as part C. A plot of the real triple product function does not appear to be significantly different than this plot. The amplitude scale here and in panel E is linear. (F) However, when the bicoherence, BIC, is plotted here, phase coherence features are visible in a wide distribution of frequency.
Clinical Applications of Frequency Domain Methods
Early in his survey of human EEG, Hans Berger identified several generic EEG patterns that were loosely correlated with psychophysiologic state. These types of activity, such as the alpha rhythms seen during awake periods with eyes closed, occurred within a stereotypic range of frequencies that came to be known as the alpha band. Eventually, five such distinct bands came to be widely accepted (Table 1).
Using an FFT it is a simple matter to divide the resulting power spectrum from an epoch of EEG into these band segments, then summate all power values for the individual frequencies within each band to determine the band power. Relative band power is simply band power divided by power over the entire frequency spectrum in the epoch of interest.
In the realm of anesthesia-related applications, traditional band power analysis is of limited utility because the bands were defined for the activity of the awake or natural sleep-related EEG without regard for the altered nature of activity during anesthesia. Drug-related EEG activity can often be observed to pass smoothly between bands as the dose changes. Familiarity with band analysis is still necessary because of the extensive literature using it.
In an effort to improve the stability of band-related changes, Volgyesi introduced the augmented delta quotient (ADQ). [53]This value is approximately the ratio of power in the band 0.5-3.0 Hz to the power in the 0.5-30.0 Hz range. This definition is an approximation because the author used analog bandpass filters with unspecified but gentle roll-off characteristics that allowed them to pass frequencies outside the specified band limits with relatively little attenuation.
John [54]applied a normalizing transformation [38]to render the probability distribution of power estimates of the delta frequency range close to a normal distribution in the CIMON EEG analysis system (Cadwell Laboratories, Kennewick, WA). After recording a baseline “self-norm” period of EEG, increases in delta band power, which are larger than three standard deviations from the self norm, were considered to represent an ischemic EEG change. [55]Other investigators have concluded this indicator may be non-specific. [56]
Another approach to simplifying the results of a power spectral analysis is to find a parameter that describes a characteristic of the spectrum. The first of these descriptors was the peak power frequency (PPF), which is simply the frequency in a spectrum at which occurs the highest power in that epoch. The median power frequency (MPF) is that frequency which bisects the spectrum: half the power is above, half below. The spectral edge frequency (SEF) is the highest frequency in the EEG, i.e., the high frequency edge of the spectral distribution. The original SEF algorithm used a form of pattern detection on the power spectrum to emulate mechanically visual recognition of the “edge.” Beginning at 32 Hz, the power spectrum is scanned downward to detect the highest frequency where four sequential spectral frequencies all contain above a predetermined threshold of power. This approach provides more noise immunity than the alternative computation, SEF95 (unpublished results, Rampil IJ, Sasse FJ, 1977). SEF95 is the frequency below which 95% of the power in the spectrum resides. Clearly either approach to SEF calculation provides a monitor that is only sensitive to changes in the width of the spectral distribution (there is always some energy in the low-frequency range). Many commonly used general anesthetics produce burst suppression EEG patterns without slowing the waves present during the remaining bursts; thus the SEF of the epoch would not reflect the anesthetic-induced depression. Combining the SEF with the BSR parameter to form the burst-compensated SEF (BcSEF;Equation 9) creates a parameter that appears to smoothly track changes in the EEG because of either slowing or suppression from isoflurane or desflurane. [11,48]Equation 9
Spectral QEEG parameters like MPF or SEF compress into a single variable the 60 or more spectral power estimates that constitute the typical EEG spectrum. As Levy pointed out, [57]a single feature may not be sensitive to all possible changes in spectral distribution. However, there is no evidence, at present, suggesting that additional parameters (describing a complex spectrum) improve the clinical utility of simple univariate parameters. Frequency domain-based QEEG parameters, like their time domain-based relatives, are generally averaged over time before display. The author uses non-linear smoothing when computing SEF that strongly filters small variations but passes large changes with little filtering. This approach decreases noise, but briskly displays major changes, such as may occur secondary to ischemia or after bolus injection of anesthetics.
The quantitative EEG variables described to this point were all created to measure patterns apparent by visual inspection in the raw waveform or the power spectrum of the EEG. Although many of these QEEG variables detect changes in the EEG caused by anesthetic drugs, all suffer from the inability to be calibrated to behavioral endpoints. Their performance as anesthetic monitors also suffers because of their sensitivity to the different EEG patterns induced by different drugs.
Bispectral Index
The recently introduced bispectral index (BIS, Aspect Medical Systems, Natick, MA) is a complex parameter, composed of a combination of time domain, frequency domain, and high order spectral subparameters (described as features of the EEG in previous literature). It is unique among the QEEG parameters reviewed here because it integrates several disparate descriptors of the EEG into a single variable based on a large volume of clinical data to synthesize a combination that correlates behavioral assessments of sedation and hypnosis yet is insensitive to the specific anesthetic or sedative agents chosen. Further, devices that implement BIS are the only ones currently approved by the Food and Drug Administration for marketing as monitors of anesthetic effect on the brain. The particular (proprietary) mixture of subparameters in BIS version 3 was derived empirically [58,59]from a prospectively collected database of EEG and behavioral scales (Figure 14), representing approximately 1,500 anesthetic administrations ([almost equal to] 5,000 h of recordings) that used a variety of anesthetic protocols. BIS was then tested prospectively in other populations. As noted, BIS has been empirically demonstrated to correlate with behavioral measures of sedation and light anesthesia [14,15,19,24,25,60]because of a variety of anesthetics, and some have suggested that close titration of anesthetic effect may improve some measures of patient outcome or operating suite efficiency. [61]
Figure 14. The development of the BIS by Aspect Medical Systems, Inc., proceeded in a stepwise fashion. First, a library of artifact-free EEG (with concurrent behavioral correlates) was accumulated. A range of prospective subparameters were calculated, and their correlation with behavior was tested. The parameters with the best performance were entered into a multivariate analysis for the creation of a final composite parameter, the BIS. The performance of the BIS has been enhanced by an iterative process that involved at least three major passes of new data collection, modeling, and progressive refinement.
Figure 14. The development of the BIS by Aspect Medical Systems, Inc., proceeded in a stepwise fashion. First, a library of artifact-free EEG (with concurrent behavioral correlates) was accumulated. A range of prospective subparameters were calculated, and their correlation with behavior was tested. The parameters with the best performance were entered into a multivariate analysis for the creation of a final composite parameter, the BIS. The performance of the BIS has been enhanced by an iterative process that involved at least three major passes of new data collection, modeling, and progressive refinement.
The calculation of BIS (Figure 15) begins with sampled EEG that is filtered to exclude high- and low-frequency artifacts and divided into epochs of 2-s duration. A series of algorithms next detect and attempt to remove or ignore artifacts. The first phase of artifact handling uses a crosscorrelation of the EEG epoch with a template pattern of an ECG waveform. If ECG or pacer spikes are detected, they are removed from the epoch, and the missing data are estimated by interpolation. Epochs repaired in this phase are still considered viable for further processing. Next eyeblink events are detected, again relying on their steroetypical shape to match a template with crosscorrelation. Epochs with blink artifact are considered to have unrepairable noise and are not processed further. Surviving epochs are checked for wandering baseline (low-frequency electrode noise), and if this state is detected, additional filtering to reject low frequencies is applied. In addition, the variance (i.e., the second central moment) of the EEG waveform for each epoch is calculated. If the variance of an epoch of raw EEG changes markedly from an average of recent previous epochs, the new epoch is marked as “noisy” and not processed further; however, the new variance is incorporated into an updated average. If the variance of new incoming epochs continues to be different from the previous baseline, the system will slowly adapt as the previous average changes to the new variance. Presuming the incoming EEG epoch is artifact-free or is deemed repaired, the time domain version of the epoch is used to calculate the degree of burst suppression with two separate algorithms: BSR and “QUAZI”(Table 2). The BSR algorithm used by the BIS calculation is similar to that described in the previous section. The QUAZI suppression index was designed to detect burst suppression in the presence of wandering baseline voltage. QUAZI incorporates slow wave (< 1.0 Hz) information derived from the frequency domain to detect burst activity superimposed on these slow waves, which would “fool” the original BSR algorithm by exceeding the voltage criteria for electrical “silence.” The waveform data in the current epoch are prepared for conversion to the frequency domain by a Blackman window function, as illustrated in Figure 10A; then the FFT and the bispectrum of the current EEG epoch are calculated. The resulting spectrum and bispectrum are smoothed using a running average against those calculated in the previous minute; then the frequency domain-based subparameters “SynchFastSlow” and “BetaRatio” are computed. The BetaRatio subparameter is the log ratio of power in two empirically derived frequency bands: log((P30-47Hz)/(P11-20Hz)). The SynchFastSlow subparameter is the contribution from bispectral analysis. SynchFastSlow is defined as another log ratio. Here the log of the ratio of the sum of all bispectrum peaks in the area from 0.5 to 47 Hz over the sum of the bispectrum in the area 40-47 Hz. The resulting BIS is defined as a proprietary combination of these QEEG subparameters. Each of the component subparameters was chosen to have a specific range of anesthetic effect where they perform best, i.e., the SynchFastSlow (HOS) parameter is well correlated with behavioral responses during moderate sedation or light anesthesia. The combination algorithm that determines BIS therefore weights the Beta Ratio (FFT) most heavily when the EEG has the characteristics of light sedation. The SynchFastSlow (bispectral component) predominates during the phenomenon of EEG activation (excitement phase) and during surgical levels of hypnosis, and the BSR and QUAZI detects deep anesthesia. The subparameters are combined using a nonlinear function whose coefficients were determined by the process described in Figure 14. Representative data illustrating the complete process of BIS determination are shown in Figure 16. Two key features of the Aspect multivariate model are that it accounts for the nonlinear stages of EEG activity by allowing different subparameters to dominate the resulting BIS as the EEG changes its character with increasing anesthesia; and second, the model framework is extensible, so new subparameters can be added to improve performance, if needed, in the presence of new anesthetic regimens. The combination of the four sub-parameters produces a single number, BIS, which decreases continuously with decreasing level of consciousness (hypnosis). As described previously, computation of bispectral parameter requires averaging several epochs; therefore, the BIS value reported on the front panel of the A-1050 monitor represent an average value derived from the previous 60 s of useable data.
Figure 15. Flowchart for the calculation of the bispectral index (BIS).
Figure 15. Flowchart for the calculation of the bispectral index (BIS).
Figure 16. Representative data from a single human volunteer demonstrating changes in EEG with increasing serum concentrations of propofol. In each of the four concentrations, 4 s of raw EEG is plotted in the top half of the figure, along with a representative time domain parameter (root mean square amplitude). In this subject, no burst suppression activity was seen, so the SR and QUAZI values are zero and not displayed. In the lower left is the corresponding power spectrum and spectral parameters computed from the same EEG epoch. The beta ratio is a component of the BIS. The lower right quadrant displays the bispectrum, again computed from the same raw data (and averaged with the previous 56 s of EEG). The SynchFastSlow component and the final BIS are also listed. These data and plots are courtesy of J. Sigl, Ph.D., Aspect Medical Systems, Inc.
Figure 16. Representative data from a single human volunteer demonstrating changes in EEG with increasing serum concentrations of propofol. In each of the four concentrations, 4 s of raw EEG is plotted in the top half of the figure, along with a representative time domain parameter (root mean square amplitude). In this subject, no burst suppression activity was seen, so the SR and QUAZI values are zero and not displayed. In the lower left is the corresponding power spectrum and spectral parameters computed from the same EEG epoch. The beta ratio is a component of the BIS. The lower right quadrant displays the bispectrum, again computed from the same raw data (and averaged with the previous 56 s of EEG). The SynchFastSlow component and the final BIS are also listed. These data and plots are courtesy of J. Sigl, Ph.D., Aspect Medical Systems, Inc.
Figure 17. Figure 16(continued).
Divining the anesthetic effect message within the EEG has long been sought. The recent availability of a device that claims to do just that is a step forward in this endeavor. Although there is currently no theoretical or mechanistic link proposed between neural network physiology in the cerebral cortex and the intrafrequency coupling notion of the BIS, the empirical correlation appears to exist. The exact role and limitations of this new technology will be determined by additional clinical experience. With the attainment of this present benchmark level of clinical correlation, further refinements in signal processing can now be reasonably expected to create increasingly useful tools for a range of clinical settings.
The authors thanks J. Sigl, Ph.D., of Aspect Medical Systems, Inc., for his thoughtful comments on the manuscript and for providing Figure 16.
[dagger]Even the 50 or 60 Hz powerline signal has potential value: it may be analyzed to assess the quality of the electrode connections.