PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Neural Eng. Author manuscript; available in PMC Oct 9, 2013.
Published in final edited form as:
PMCID: PMC3793904
NIHMSID: NIHMS523716
Burst suppression probability algorithms: state-space methods for tracking EEG burst suppression
Jessica Chemali,1 ShiNung Ching,1,2 Patrick L Purdon,1,2 Ken Solt,1,2 and Emery N Brown1,2,3
1Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital, Boston, MA, USA
2Department of Brain and Cognitive Science, Massachusetts Institute of Technology, Cambridge, MA, USA
3Institute for Medical Engineering and Science, Massachusetts Institute of Technology, Cambridge, MA, USA
Emery N Brown: enb/at/neurostat.mit.edu
Objective
Burst suppression is an electroencephalogram pattern in which bursts of electrical activity alternate with an isoelectric state. This pattern is commonly seen in states of severely reduced brain activity such as profound general anesthesia, anoxic brain injuries, hypothermia and certain developmental disorders. Devising accurate, reliable ways to quantify burst suppression is an important clinical and research problem. Although thresholding and segmentation algorithms readily identify burst suppression periods, analysis algorithms require long intervals of data to characterize burst suppression at a given time and provide no framework for statistical inference.
Approach
We introduce the concept of the burst suppression probability (BSP) to define the brain’s instantaneous propensity of being in the suppressed state. To conduct dynamic analyses of burst suppression we propose a state-space model in which the observation process is a binomial model and the state equation is a Gaussian random walk. We estimate the model using an approximate expectation maximization algorithm and illustrate its application in the analysis of rodent burst suppression recordings under general anesthesia and a patient during induction of controlled hypothermia.
Main result
The BSP algorithms track burst suppression on a second-to-second time scale, and make possible formal statistical comparisons of burst suppression at different times.
Significance
The state-space approach suggests a principled and informative way to analyze burst suppression that can be used to monitor, and eventually to control, the brain states of patients in the operating room and in the intensive care unit.
Burst suppression is an electroencephalogram (EEG) pattern indicating a severe reduction in the brain’s neuronal activity and metabolic rate [1]. Observed in profound general anesthesia [2], coma [3], hypothermia [4], epilepsy due to Ohtahara’s syndrome [5], and postasphyctic newborns [6], burst suppression consists of periods of electrical activity alternating with a flat line or isoelectric state termed a suppression. Both the bursts and the suppression periods can last from a few seconds to several minutes [7]. Figures 1(A) and (C) show respectively a 5-min and a 1-min segment of burst suppression induced by the anesthetic propofol.
Figure 1
Figure 1
A. 5-min of burst suppression recorded from a patient following a propofol bolus. The 5-microvolt threshold used for the detection of the binary events is shown in red. B. The corresponding binary time-series where 1 represents a suppression and 0 represents (more ...)
Devising accurate, reliable methods to quantify burst suppression is an important clinical and research problem. Medical coma is often induced by administering an anesthetic such as propofol for cerebral protection following a brain injury or to arrest intractable epilepsy [4, 8]. The level of burst suppression is continuously monitored as a marker of the level of coma in order to balance the trade-offs between the drug‘s therapeutic benefits and its side effects [9]. Induced hypothermia is used for brain protection in patients recovering from a cardiac arrest and in patients having certain types of cardiac, brain and major vascular surgeries [10]. The ascent into and the descent out of hypothermia can be tracked by monitoring the change in temperature along with the change in burst suppression [11]. Persistence of burst suppression patterns in the EEG of patients in coma is commonly associated with a poor prognosis [3]. EEG-based monitors used to track the brain states of patients under general anesthesia compute measures of burst suppression as part of their analysis algorithms [12, 13].
In anesthesiology research, monitoring the level of burst suppression has been proposed as a way of assessing cortical reactivity to stimuli and hence for making inferences about the sites and mechanisms of anesthetic drug actions [14]. The effectiveness of anesthetic arousal agents has been tested by measuring their ability to induce emergence from burst suppression [15, 16]. Finally, one way of characterizing the potential utility of a new anesthetic drug is by measuring its efficacy in maintaining a specified level of burst suppression using a computer-controlled infusion with feedback [17, 18].
Quantification of burst suppression begins by thresholding and segmenting the EEG [1823]. Thresholding sets a voltage level to separate burst and suppression events. If the EEG is less than the threshold in the interval, the event is a suppression and is assigned a value of 1 whereas, if the EEG is greater than the threshold in the interval then the event is a burst and is assigned a value of 0. Figures 1(B) and (D) show the binary time-series for a threshold of 5 microvolts and for the EEG signals in figures 1(A) and (C). Using the binary time-series, burst suppression is commonly characterized by computing the burst suppression ratio (BSR); the fraction of a given time interval that the EEG is suppressed [24]. The BSR is a number ranging from 0, meaning no suppression to 1, meaning an isoelectric EEG.
The BSR is positively correlated with reduction in cerebral metabolic rate (CMR) [25]. During general anesthesia and during induced hypothermia when this fraction increases to one, the CMR decreases in a dose-dependent manner until it plateaus [9, 2628]. Other approaches to characterizing burst suppression have included entropy measures [29, 30] and artificial neural networks and support vector machines [31].
Although the importance of quantitatively analyzing burst suppression is broadly appreciated, there are key shortcomings with current approaches. The binary time-series can be computed on intervals as short as 100 msec, yet BSR estimation requires several consecutive seconds or minutes of binary events [15, 18]. This assumes that the brain state remains constant during these long estimation intervals, a condition which does not hold during the transitions into or out of general anesthesia or hypothermia. An important objective of quantitative burst suppression analyses is to make formal statistical comparisons at different time points. The statistical properties of the BSR estimated by averaging the binary events over long time intervals have not been described. Hence, it is not clear how best to use current BSR estimates in formal statistical analyses of burst suppression.
We present a new approach to conducting dynamic analyses of burst suppression based on the state-space framework for point processes and binary observations developed in [3234]. The observation model is a binomial process and the temporal evolution of the brain state of burst suppression is defined by a state equation represented as a Gaussian random walk. By making a logistic transformation on the state, we introduce the concept of the burst suppression probability (BSP) to define the brain’s instantaneous probability of being in the suppressed state. We estimate the model parameters using an approximate EM algorithm and illustrate its application in a rodent study of the effects of the cholimimetic drug physostigmine on burst suppression, and in a study of burst suppression induced in a patient by controlled hypothermia. Our approach obviates the need to average binary events over long intervals and allows formal statistical comparisons of burst suppression at different time points.
2.1. Model definition
To formulate our state-space model, we follow the state-space paradigm for analyzing point processes and binary time-series in [3234]. We assume that the EEG recordings segmented into binary events are collected in the observation interval (0, T] and that our state-space model is defined on a discrete set of lattice points within that interval. To define the lattice, we choose I large and divide (0, T] into I subintervals of equal width Δ = T I−1. The state-space model is evaluated at iΔ for i = 1…I.
A state-space model is characterized by its state and observation equations. The state equation defines the unobservable state process whose evolution we wish to track over time. In our case, the state represents the brain’s state of burst suppression. We define our state to be positively related to the probability of suppression. That is, as the state increases the probability of suppression increases and as the state decreases the probability of suppression decreases. The observation equation describes how the observations relate to the unobservable state process. Our objective is to estimate the brain’s burst suppression state, its instantaneous probability of being suppressed and the associated confidence intervals.
We assume that in each interval Δ is divided into n subintervals of width δ and so that Δ = and in each subinterval there can be at most one suppression event or one burst event. Let bi be the number of suppression events in iΔ. We assume that the observation model is described by the binomial probability mass function as
equation M1
(1)
where pi, the BSP, is
equation M2
(2)
where xi is the brain’s burst suppression state at time iΔ. In other words, pi is the instantaneous probability of being suppressed. The logistic function (2) links the brain’s burst suppression state to the probability of a suppression event and ensures that pi remains between 0 and 1 as xi ranges across all real numbers.
We define the state equation as the random walk
equation M3
(3)
where the εi are independent, Gaussian random variables with mean 0 and variance equation M4. This definition of the state provides a stochastic continuity constraint which ensures that the states, and hence that the BSPs that are close in time are close in value. The parameter equation M5 governs how rapidly the BSP can change; the larger (smaller) the value of equation M6 the more rapidly (slowly) the state and the BSP can change.
2.2. Model estimation
To present our estimation algorithm we take b = (b1, b2, …, bI) and x = (x1, x2, …, xI). Based on the random walk defined in (3), the joint probability density of the state process is
equation M7
(4)
given the initial state x0. The joint probability mass function of the observed suppression events given the states is
equation M8
(5)
Our objective is to estimate using maximum likelihood (ML) the state process x and the parameters equation M9 and x0, where we treat x0 as a parameter. Once we obtain these estimates we can readily compute the BSP with its confidence intervals.
To compute the ML estimates of the parameters we use an approximate expectation maximization (EM) algorithm for point processes and binary time-series [3234]. The EM algorithm is a well established method for simultaneously estimating model parameters and an unobservable state process by iteratively maximizing the expectation of the complete data log likelihood [35]. The complete data likelihood is
equation M10
(6)
The EM algorithm is
Expectation step
At iteration [ell] + 1 we compute in the expectation step the expected value of the complete data log likelihood given the data b and the estimates equation M11 and equation M12 of the parameters from iteration [ell]:
equation M13
(7)
We see upon expanding the right side of (7) that we need to estimate three quantities for i = 1,…, I. The expectations of the state variables conditioned on the data up to time I are
equation M14
(8)
and the covariances of the state variables conditioned on the data up to time I are
equation M15
(9)
and
equation M16
(10)
To compute these quantities efficiently, we divide the expectation step into three parts. First we compute the estimates of xi|i and equation M17 using a binary filter algorithm. Second, we use a fixed interval smoothing (FIS) algorithm to compute xi|I and equation M18. Finally, we use the state-space covariance algorithm to compute the covariances Wi|I and Wi,i−1|I.
Binary filter algorithm
Given the parameter estimates from iteration [ell], this step estimates xi|i and equation M19, the state and the variance at i, given the data from the start of the observation interval up through time i [32, 36].
The one step prediction mean and variance are given by
equation M20
(11)
equation M21
(12)
The posterior mode and variance are given by
equation M22
(13)
equation M23
(14)
The initial conditions are equation M24 and pi|i is (2) evaluated at xi|i. This filter is nonlinear because xi|i appears on both sides of (13). We can compute it using Newton’s method [32, 36] However, when Δ is small, adjacent states are close, and we can replace the term pi|i in (13) with pi−1|i−1.
Fixed interval smoothing algorithm
Using the estimates from the binary filter, a FIS algorithm [32, 33] gives the state and variance estimates xi|I and equation M25 respectively for i = I − 1, …, 1. They are the estimates at time i conditioned on all the data up through time I. The final state estimate is thus a Gaussian variable with mean xi|I and variance equation M26. The FIS is
equation M27
(15)
equation M28
(16)
equation M29
(17)
The initial conditions are xI|I and equation M30, computed at the final step of the binary filter algorithm.
State-space covariance algorithm
Finally, using the state-space covariance algorithm [37] we compute σi, j|I :
equation M31
(18)
where 1 ≤ ijI. The covariances we require in (9) and (10) are thus given by
equation M32
(19)
and
equation M33
(20)
Maximization step
To carry out the maximization step at iteration [ell] + 1 we let equation M34 and assume that τ has a gamma prior density defined as
equation M35
(21)
for α > 1 and β > 0. We maximize the expected value of the complete data log likelihood with respect to τ using the gamma prior density for τ in (21). The expected value of the complete data log likelihood serves as the likelihood in the expression for the posterior. The log posterior density of τ is proportional to
equation M36
(22)
We chose the gamma distribution for the prior because it is flexible and because it is the conjugate distribution for τ in the Gaussian likelihood [32]. This latter property allows us to compute the update of τ in closed form in the M-step (23).
We maximize (22) with respect to τ and find
equation M37
(23)
We maximize (22) with respect to x0 and find
equation M38
(24)
The algorithm iterates between the expectation and maximization steps until convergence. The convergence criteria we use are the same as those developed by [35]. The ML estimates of τ, or equivalently equation M39, and x0 are respectively equation M40 and equation M41, where L is the last iteration of the algorithm.
When (11–14) and (15–17) are applied to data to compute pi|i and pi|I with x0 and equation M42 evaluated at their ML estimates, we term them respectively the BSP filter and smoothing algorithms.
2.3. The probability density of the BSP and confidence intervals
The FIS estimate, xi|I together with (2) gives us the probability of a suppression at time i for i = 1, …, I. Through a change of variable, we can then compute the probability density function for the BSP at time i as
equation M43
(25)
The 95% confidence intervals are computed from the cumulative density of (25) by identifying its 2.5th and 97.5th percentiles.
2.4. Comparison of BSPs at different times
Given the ML estimates of x0 and equation M44, it follows that we can compute with the binary filter, the FIS algorithm and the state-space covariance algorithm the Gaussian approximation to the joint posterior density of the states x. This provides an empirical Bayes estimate of the joint posterior density of the burst suppression states [35]. Because the logistic transformation (2) that relates the state xi to BSP pi is monotonic, we can compute the joint posterior density of p from the joint posterior density of x. Therefore, we can make formal inferences comparing the BSP at any time i, with the BSP at any time j by computing an empirical Bayes estimate of the posterior probability that pj > pi, which is equivalent to the posterior probability that xj > xi. We do so using a Monte Carlo approach [34]. Using the covariance algorithm, for times i and j such that 1 ≤ i < j the covariance between states at time i and time j is
equation M45
(26)
We can then draw M samples from the Gaussian distribution with mean
equation M46
(27)
and covariance matrix
equation M47
(28)
and count the number G of instances in which xj|I > xi|I. The estimate of the probability of interest is
equation M48
(29)
In our analysis we chose M=10 000.
3.1. Detection of the binary events
We applied a two-step procedure consisting of bandpass filtering followed by thresholding to convert the EEG recordings into time-series of binary events. For all of the experiments the EEG signals were bandpass filtered between 5 and 30 Hz. For the rat experiment the filtered signal was thresholded at 50 microvolts and segmented at 1 Hz (δ = 1s), whereas for the human experiment the threshold was set at 5 microvolts and the segmentation rate was again 1 Hz (δ = 1s). If the EEG exceeded the threshold in Δ, then the event was classified as a burst, or a 0, whereas if the EEG was less than the threshold in Δ then the event was classified as a suppression, or a 1.
3.2. The binomial observation models
For both the rat and human experiments, we analyzed the data in (Δ = 1) s intervals, making the observation model for both the rat and the human experiments the binomial with n = 1. For filtering and smoothing binary processes the update interval can be chosen arbitrarily small as is appropriate for the problem being studied as long as the state and observation processes are defined in continuous time. We chose (Δ = 1) s because in the operating room and in the intensive care unit using 1 s updates to track burst suppression balances the trade-off between unnecessarily frequent updates and missing important changes in the process due to infrequent sampling. This choice of Δ also illustrates that in the state-space framework, unlike with the BSR, large intervals are not needed to produce smooth BSP estimates.
3.3. Choice of the prior density for τ
We used an empirical Bayes’ approach to selecting the prior for τ by using the location of the likelihood to help guide the choice of the prior [35]. We first estimated τ without a prior (effectively an uninformative prior) to gain insight into its relative size and scale. Given the initial estimate, we then constructed the following prior for our analyses: for the physostigmine analysis, we chose α = 105 and β = 2 and for the hypothermia analysis we chose α = 3 × 105 and β = 2.
3.4. The BSR smoothing and filter algorithms
Based on previous reports using the BSR [17, 18] we applied 4 different types of symmetric and one-sided BSR filters for comparison with our BSP algorithms making the observation model for both the rat and the human experiments the binomial with n = 1. The symmetric BSR filters are: 15-s symmetric filter with no overlap; 15-s symmetric filter with 14-s overlap; 60-s symmetric filter with no overlap; and 60-s symmetric filter with 59-s overlap. We also used the same 15-s and 60-s bandwidths with the same degree of overlap to construct one-sided BSR filters. We define BSR symmetric and one-sided filters at time i as
equation M49
(30)
equation M50
(31)
We used the symmetric BSR filters to compare with our BSP smoothing algorithm in the analysis of the physostigmine experiment and we used the one-sided BSR filters to compare with our BSP filter algorithm in the analysis of the hypothermia experiment.
We illustrate the application of our BSP algorithms through comparisons with the BSR in two analyses: a rat in general anesthesia-induced burst suppression administered physostigmine to elicit arousal; and a patient undergoing controlled hypothermia for cerebral protection during total circulatory arrest.
4.1. Burst suppression and the arousal effects of physostigmine
Physostigmine is a cholimimetic drug that has been used in anesthesiology research to induce emergence from general anesthesia [38] and in anesthesiology practice to treat emergence delirium [39], a confusional state that some patients, frequently children, enter on emergence from general anesthesia. For both induction of emergence and treatment of emergence delirium, the effect of physostigmine is attributed to an increase in cortical levels of the excitatory neurotransmitter acetylcholine [39]. Because administration of physostigmine induces cholinergically-mediated arousal the administration of physostigmine to an animal maintained in burst suppression should induce at least a decrease in the level of burst suppression. Therefore, in this experiment, to quantify the time course of the effect of physostigmine on burst suppression we used the BSP smoothing algorithm (15–17) and the BSR symmetric filters: 15-s symmetric filter with no overlap (figure 2(C), green curve); 15-s symmetric filter with 14-s of overlap (figure 2(C), red curve); 60-s symmetric filter with no overlap (figure 2(D), green curve); and 60-s symmetric filter with 59-s of overlap (figure 2(D), red curve).
Figure 2
Figure 2
A. The EEG of an isoflurane anesthetized rat administered physostigmine to assess its arousal effects. At minute 10 (red vertical arrow) normal saline is injected. At approximately minute 16 (red star), physostigmine is injected and the EEG promptly switches (more ...)
This study was approved by the Massachusetts General Hospital Subcommittee on Research Animal Care. A rat implanted with extradural EEG electrodes was anesthetized with 2% isoflurane to induce burst suppression. This concentration was maintained for 70-min. We analyzed the EEG recorded during the last 40-min. An intravenous injection of saline was administered as a control stimulus at minute 10 (figure 2(A), vertical arrow). At minute 16, physostigmine was administered intravenously (figure 2(A), star).
Burst suppression was readily visible in the raw EEG (figure 2(A)) and in the binary time-series (figure 2(B)). As expected, injection of normal saline at minute 10 (figure 2(A), red arrow) had no effect on the raw EEG or on the binary time-series. In contrast, the effect of injecting physostigmine at minute 16 (figure 2(A), red star) was clearly evident in both series. To quantify the structure in the EEG data, we fit the BSP model to the binary series. The model fitting required 63-s. From both the BSR and the BSP estimates (figures 2(C) and (D)) it is easy to see that both the BSR and BSP (figures 2(C) and (D)) are round 0.5 with the rat receiving 2% isoflurane at the start of the experiment and that it returned to approximately this value as the effect of the physostigmine subsided.
Following the physostigmine injection at minute 16, the BSP (figure 2(D)) decreased to zero and remained at zero until minute 24. It then returned to around 0.5 at minute 31 where it stayed for the balance of the experiment. The BSP produced smooth time-series estimates of the instantaneous probabilities of burst suppression with confidence intervals (figure 3(A)). Although the inspired concentration of isoflurane was maintained at 2%, the BSP fluctuated between 0.4 and 0.6 suggesting that the pharmacokinetic state defined by a fixed inspired concentration does not agree necessarily with the neurophysiological state of the brain.
Figure 3
Figure 3
A. BSP smoothing algorithm estimate (black curve) and its associated 95% confidence (Bayesian credibility) intervals (red curves) for minutes 2–5 in figure 2(C). B. BSR estimate computed with a 60-s interval and 59-s overlap (black curve) and (more ...)
The BSR computed from the 15-s windows without (figure 2(C), green) and with (figure 2(C), red) 14-s overlap showed much more variability than the BSP. The overlapping filter estimates added high frequency noise to the BSR smoothing procedure. The non-overlapping BSR estimate is a sub-sample of the overlap BSR estimate so as expected, it was less noisy and the two estimates agreed every 15-s. In contrast, the BSR estimate computed from the 60-s windows without overlap (figure 2(D), green) and with (figure 2(C), red) 59-s overlap showed much less variability and agreed more closely with the time course of the BSP. The non-overlapping 60-s smoother oversmoothed the data as features that are readily apparent in the overlapping smoother and the BSP, such as the fluctuations between minute 26 and minute 29 and between minute 31 and minute 36, were lost with the non-overlapping smoother (figure 2(D)). The BSR estimated from the 15-s windows suggested that during the experiment with a constant inspired concentration of isoflurane, the probability of burst suppression fluctuated between 0.15 and 0.95, whereas the degree of fluctuation estimated by the 60-s window BSR estimates agreed more closely with that computed from the BSP.
Although the BSP stayed around 0.5 before and after the physostigmine injection, it shows fluctuations between 0.4 and 0.6 (figures 2(C) and 3(A)). For each BSP estimate we can compute approximate 95% confidence intervals based on (25). The benefit of the 95% confidence intervals is that they give a measure of the uncertainty in the BSP estimate and they allow us to infer whether the BSP at two pre-selected times, say at minute 2.5, where the BSP is 0.55, and at minute 4.5, where the BSP is 0.64 differ (figure 3(A)). In this case, the estimates at these two points did differ significantly as suggested by the fact that the respective 95% confidence intervals did not overlap. This is verified by the fact that the 95% confidence interval for the difference (0.09 ± 0.04 = [0.05, 0.13]) does not include zero.
The BSR estimates at minute 2.5 and minute 4.5 are 0.49 and 0.60, respectively. For comparison, we computed, using the Gaussian approximation to the binomial [31], 95% confidence intervals based on the BSR estimates (figure 3(B)). These intervals are narrower than the ones derived from the BSP. We also computed the 95% confidence interval for the difference between the two BSR estimates at minute 2.5 and 4.5 (0.11 ± 0.008 = [0.082, 0.098]) which did not include zero and was narrower than the 95% confidence interval for the difference based on the BSP. The BSR confidence interval underestimates the uncertainty in the analysis because it assumes that the observations are independent, whereas the BSP algorithm models the dependence in the binary time-series through the state-space model.
Using our methods to make statistical inferences about the effect of physostigmine on burst suppression, not only at pre-selected time points but across the entire experiment, is an important feature of our framework. Our BSP algorithm estimates the joint distribution of the state process thus, we can evaluate the empirical Bayes, or equivalently ML, estimate of Pr(xj > xi) = 0.975 for any 0 ≤ i < jI. This is equivalent to the probability that the BSP at time j is greater than the BSP at time i because the transformation between the state variable xj and the BSP pj is monotonic (2). We can therefore make formal comparisons among the BSPs and state when any BSP value differs from another.
These I(I − 1)/2 comparisons are easily represented in a lower triangular matrix in which every time point along the horizontal axis is compared to every preceding time point (figure 3(C)). A red entry corresponds to Pr(pj > pi) > 0.975, a black entry corresponds to Pr(pj > pi) < 0.025 and a gray entry corresponds to 0.025 ≤ Pr(pj > pi) ≤ 0.975 where i is the index across the horizontal axis and j is the index across the vertical axis.
The change in the BSP following the physostigmine injection was dramatic, and could be easily seen in the raw data (figure 2(A)). However, the point-by-point comparison matrix confirms that following injection of physostigmine at minute 16 until about minute 26 (figure 3(C), red area, minutes 16–26, x-axis) the BSP was significantly smaller than at the time points prior to the injection. The coordinate 25, 20 in figure 3(C) is black because P25 and P20 are both close to 0. The gray patterns pre- and post-saline injection (figure 3(C), minute 10, x-axis) are similar suggesting, as expected, that saline had no effect. These results show that our BSP state-space model provides a principled statistical approach to measuring quantitatively the effect of physostigmine on BSP. These analyses are not possible with the BSR as these algorithms do not consider the joint distribution of their estimates across time.
4.2. Burst suppression and hypothermia
As stated in the Introduction, burst suppression can be produced by states of profound general anesthesia and hypothermia. Hypothermia is commonly induced in patients having cardiac and major vascular surgery for cerebral protection [9]. To illustrate the relevance of our algorithms for real-time monitoring we analyzed the EEG recordings of a patient under general anesthesia who underwent controlled hypothermia as part of total circulatory arrest for thoracic aortic aneurysm repair.
This study was approved by the Massachusetts General Hospital Human Research Committee. Written informed consent was not required as the EEG data are standard measurements recorded as part of standard anesthesia care in our institution. We fit the state-space model to the patient’s binary time-series derived from the scalp EEG recorded from a Sedline (Masimo, Irvine, CA) monitor with its standard six-electrode montage during the transition into hypothermia (figure 4(A)). The entire data set consisted of 208-min of EEG recordings. We analyzed the patient’s 45-min transition into the isoelectric state (figure 4(A)) using the BSP filter (figures 4(C) and (D), black curve). The patient was in a standard state of general anesthesia appropriate for surgery from minute 0 until minute 8. As part of the procedure to induce total circulatory arrest, hypothermia was initiated at minute 8 to cool the patient to 19 °C. During approximately the next 25 min, the EEG evolved from that observed in a standard surgical state of general anesthesia through burst suppression (figure 4(A), minutes 8–37). Between minutes 37–45, the EEG was isoelectric. These transitions observed in the raw EEG were also evident in the binary time-series (figure 4(B)).
Figure 4
Figure 4
A. The EEG recorded in a patient undergoing controlled hypothermia. B. The binary time-series associated with A. C. The BSP filter estimate (black curve), the one-sided BSR estimate computed using 15-s intervals with no overlap (green curve), and the (more ...)
To emulate real-time analysis, we analyzed these data with the BSP filter algorithm (11–14) (figures 4(C) and (D), black curve) and with 4 one-sided BSR filters: 15-s one-sided filter with no overlap (figure 4(C), green curve); 15-s one-sided filter with 14-s overlap (figure 4(C), red curve); 60-s one-sided filter with no overlap (figure 4(D), green curve); and 60-s one-sided filter with 59-s overlap (figure 4(D), red curve). We fit the BSP model to the binary observations in 64-s. The BSP filter (figures 4(C) and (D), black curve) showed quantitatively that as the temperature dropped, the EEG transitioned gradually into a deeper state of burst suppression. Between minutes 10 and 15, the BSP increased almost monotonically from 0 to approximately 0.8, indicating that the patient rapidly reached a deep level of burst suppression. The remaining transition between minutes 15 and 37 to an isoelectric state occurred with a series of what resemble logarithmic increases with rapid decreases that became progressively smaller. The dynamics in the time course of the BSP agreed closely with the patterns in the raw EEG (figure 4(A)) and the binary time-series (figure 4(B)).
Both 15-s BSR filters oscillated almost between 0 and 1 continually during the transition into the isoelectric state (figure 4(C)), giving the impression that the patient’s brain’s electrical activity moved back and forth across the entire dynamic range of burst suppression. The 60-s BSR filters agreed closely with the BSP between minute 10 and minute 14 (figure 4(D)). Between minute 14 and minute 35, these BSR filters showed an oscillatory pattern that resembled the BSP but with wider excursions. As in the previous example, the excursions of the overlapping 60-s BSR filter were greater than those of the non-overlapping 60-s filter. The overlapping filter (figure 4(D), red curve) estimated that the patient was intermittently in an isoelectric state (i.e. BSP = 1) for 15–30 s from minute 17 on. The non-overlapping filter estimated that the patient was in an isoelectric state between minutes 16 and 21. Following this point, this filter tracked the BSP. The 60-s BSR filter estimates agreed with the BSP estimates more closely than the 15-s BSR filters. However, the BSP filter provided the more informative characterization of the dynamics of the patient’s burst suppression during the transition into an isoelectric state induced by hypothermia.
Our analysis demonstrates that given estimates of the state-space model parameters, the BSP filter can be used in real-time to track burst suppression dynamics with a 1-s resolution.
Burst suppression is a state of profound brain inactivation that appears in several drug-induced and pathological conditions. We have formulated the problem of analyzing burst suppression as a dynamic signal processing question and presented a state-space model to characterize its temporal evolution. The observation model is a binomial process (1) and the state equation is a Gaussian random walk model (3). We introduced the concept of the BSP (2) as a principled way to define the instantaneous probability of the EEG being suppressed.
We estimated the state and model parameters by modifying the approximate EM algorithm for state-space estimation for binary and point processes developed in [32] to include a gamma prior distribution on the inverse of the state variance (21–22). Our approach allows us to estimate the BSP on a second-to-second time scale and to make formal statistical comparisons of burst suppression activity at different time points by computing confidence intervals and/or empirical Bayes posterior probabilities. We illustrated our new BSP algorithms in comparison to the BSR algorithms in one experimental and one clinical application.
Our state-space model approach offers several advantages over current methods for analyzing burst suppression. First, the state-space model provides a clear definition of the BSP as the instantaneous probability of being suppressed (2). Second, the BSR does not have a principled method for selecting the bandwidth and degree of overlap for its filters. We used 15-s and 60-s windows to compute the BSR because these non-overlapping filters had been used in previous reports [15, 18]. In addition, we computed the BSR estimates with overlap to provide BSR updates that matched the 1-s updates we computed from our BSP algorithms. We further constructed both symmetric and one-sided versions of both BSR algorithms to compare directly with our BSP smoothing and filter algorithms respectively. Although BSR estimates computed in 4-s intervals have also been reported [24], we did not show analyses with those estimates because they did not differ appreciably from the binary time-series. As we demonstrated, it is possible to compute smoother (rougher) estimates of the BSR by taking longer (shorter) computation windows and/or by not allowing (allowing) adjacent windows to overlap. A bandwidth selection procedure may offer one solution to guide these decisions [40].
Third, our state-space framework addresses these issues by using the state-model to impose a temporal continuity constraint on the relation of the BSP values at nearby time points. The state-space variance equation M51 governs the degree of smoothing in the BSP estimates. Larger (smaller) values of equation M52 allow for less (more) smoothness in the BSP time course. Placing a prior distribution on equation M53 places a constraint on the degree of smoothness that can be imposed by this parameter.
Although our BSP algorithm uses an empirical Bayes’ procedure to choose equation M54, the algorithm’s local prediction-and-correction scheme is another important feature that helps explain its good performance. Equation (13) shows that the update xi|i is computed based on the previous update xi−1|i−1 so that when the update interval is small, xi−1|i−1 gives a good guess of where the next state estimate is likely to be. This is the algorithm’s prediction term. The binomial innovations term, binpi|i, is the difference between the number of suppression events that is observed and the number that would be expected in the current observation interval based on the current estimate of pi|i. This term is bounded between −n and n. The left extreme occurs if the BSP is close to one and no suppression is observed, whereas the right extreme occurs if the BSP is close to 0 and n suppressions are observed. These rare events provide the maximum possible innovation or local correction to the new state estimate. The closer (more distant) the observed number of suppression events is from the prediction, the less (greater) the correction that is made to xi−1|i−1 to compute xi|i.
The term equation M55, which is the gain in the BSP filter, governs how much the innovation is weighted in computing the new update. Because equation M56 is the one-step prediction variance, the greater (less) this variance is, the greater (less) the innovation is weighted. Under our Gaussian approximation to the state, the FIS algorithm (15–17) provides an approximately optimal strategy for computing from the filter estimates state estimates that depend on all of the binary observations [32]. These local adaptive features of the BSP algorithms, which are characteristic of Kalman filter-like algorithms [41], are another reason that these BSP algorithms could be expected to perform better than the BSR methods that use only elementary filtering strategies. We have previously demonstrated that our binary smoothing algorithm performs better than ad hoc smoothing methods [33], and as well as or better than more elaborate smoothing algorithms that have an automatic bandwidth selection criterion [42].
Fourth, a significant benefit of our framework is the ability to use the model to make statistical inferences about the character of burst suppression being studied. This is especially important in studies such as the rat example in which key questions are measuring the second-to-second arousal effect of physostigmine on burst suppression and comparing the level of burst suppression before and after drug administration. Our state-space modeling framework provides an empirical Bayes’ estimate of the joint posterior distribution of the BSP estimates across time. Using Monte Carlo methods we can easily compute confidence statements (figure 3(A)) and posterior probabilities (figure 3(C)) for any functions of interest. For example, if in the physostigmine experiment we wanted to compare BSP in a pre-treatment interval with the BSP in a post-treatment interval, we can use the Monte Carlo algorithm to make pairwise comparisons of points chosen at random from the two intervals. The posterior probability that the BSP on the pre-treatment interval is greater than the BSP on the post-treatment interval is the fraction on the pairwise comparisons in which the pre-treatment point exceeded the post-treatment point. Because our inferences are based on an estimate of the joint posterior distribution of the state variables, we obviate the problems of multiple comparisons that are common to hypothesis-testing approaches in multivariate analyses.
We demonstrated how local 95% confidence intervals can be computed from the BSR estimates using the well-known Gaussian approximation to the binomial [32]. Such intervals have not been previously reported. Because these confidence intervals, unlike those computed from our BSP algorithm, are local they do not use all of the data. Moreover, because they assume that the observations are independent the BSR confidence intervals understate the uncertainty in the BSR estimates. In contrast, the BSP algorithm reports wider confidence intervals because the state-space formulation models the temporal dependence in the binary time-series.
Finally, our last example shows that given estimates of the model parameters, the BSP filter algorithm could be combined with a thresholding and segmenting algorithm to track burst suppression in patients in real-time. In this situation the model parameter estimation would be conducted either off line or on a slower time scale than that for updating the BSP estimates. Possible applications of such a real-time BSP algorithm include tracking burst suppression during surgery [13] as well as monitoring the state of medically-induced coma in patients in the intensive care unit [8, 9]. We were unable to compare our BSP filter algorithm directly with the BSR algorithms used in current depth-of-anesthesia monitors because they are proprietary [13]. Our results suggest that our algorithm should compete favorably with these procedures.
In summary our state-space paradigm for the analysis of burst suppression could be applied in research and clinical analyses of this important brain state. Studies using our paradigm to track burst suppression in real-time in the operating room and in the ICU, to study the efficacy of physostigmine and other stimulants in inducing emergence from burst suppression, to track the state of the brain in postasphytic neonates, to track the state of burst suppression in patients receiving anesthetics for maintenance of medical coma [4345] will be the topics of future reports.
Acknowledgments
This work was supported by: NIH Awards DP1-OD003646 (ENB), DP2-OD006454 (PP), K08-GM094394 (KS), and the Burroughs-Wellcome Fund Award 1010625 (SC).
Footnotes
Some figures may appear in colour only in the online journal
1. Amzica F. Basic physiology of burst-suppression. Epilepsia. 2010;50:38. [PubMed]
2. Brown EN, Lydic R, Schiff ND. General anesthesia, sleep, and coma. N Engl J Med. 2010;363:2638–50. [PMC free article] [PubMed]
3. Young GB. The EEG in coma. J Clin Neurophysiol. 2000;17:473. [PubMed]
4. Marion DW, Penrod LE, Kelsey SF, Obrist WD, Kochanek PM, Palmer AM, Wisniewski SR, DeKosky ST. Treatment of traumatic brain injury with moderate hypothermia. N Engl J Med. 1997;336:540–6. [PubMed]
5. Ohtahara S, Yamatogi Y. Epileptic encephalopathies in early infancy with suppression-burst. J Clin Neurophysiol. 2003;20:398. [PubMed]
6. Löfhede J. The EEG of the Neonatal Brain: Classification of Background Activity. Goteborg: Chalmers University of Technology; 2009.
7. Van De Velde M. Signal Validation in Electroencephalography Research. Eindhoven: Technische Universiteit Eindhoven; 2000.
8. Rossetti AO, Reichhart MD, Schaller MD, Despland PA, Bogousslavsky J. Propofol treatment of refractory status epilepticus: a study of 31 episodes. Epilepsia. 2004;45:757–63. [PubMed]
9. Doyle PW, Matta BF. Burst suppression or isoelectric encephalogram for cerebral protection: evidence from metabolic suppression studies. Br J Anaesth. 1999;83:580–4. [PubMed]
10. Stecker MM. Neurophysiology of surgical procedures for repair of the aortic arch. J Clin Neurophysiol. 2007;24:310–5. [PubMed]
11. Stecker MM, Cheung AT, Pochettino A, Kent GP, Patterson T, Weiss SJ, Bavaria JE. Deep hypothermic circulatory arrest: II. Changes in electroencephalogram and evoked potentials during rewarming. Ann Thorac Surg. 2001;71:22–8. [PubMed]
12. Soehle M, Ellerkmann RK, Grube M, Kuech M, Wirz S, Hoeft A, Bruhn J. Comparison between bispectral index and patient state index as measures of the electroencephalographic effects of sevoflurane. Anesthesiology. 2008;109:799. [PubMed]
13. Bruhn J, Bouillon TW, Shafer SL. Bispectral index (BIS) and burst suppression: revealing a part of the BIS algorithm. J Clin Monit Comput. 2000;16:593–6. [PubMed]
14. Hartikainen KM, Rorarius M, Peräkylä JJ, Laippala PJ, Jantti V. Cortical reactivity during isoflurane burst-suppression anesthesia. Anesth Analg. 1995;81:1223–8. [PubMed]
15. Luo T, Leung LS. Basal forebrain histaminergic transmission modulates electroencephalographic activity and emergence from isoflurane anesthesia. Anesthesiology. 2009;111:725. [PubMed]
16. Zhang LN, Li ZJ, Tong L, Guo C, Niu JY, Hou WG, Dong HL. Orexin-a facilitates emergence from propofol anesthesia in the rat. Anesth Analg. 2012;115:789–96. [PubMed]
17. Cotten JF, Le Ge R, Banacos N, Pejo E, Husain SS, Williams JH, Raines DE. Closed-loop continuous infusions of etomidate and etomidate analogs in rats: a comparative study of dosing and the impact on adrenocortical function. Anesthesiology. 2011;115:764. [PMC free article] [PubMed]
18. Vijn PC, Sneyd JR. IV anaesthesia and EEG burst suppression in rats: bolus injections and closed-loop infusions. Br J Anaesth. 1998;81:415–21. [PubMed]
19. Rampil IJ, et al. No correlation between quantitative electroencephalographic measurements and movement response to noxious stimuli during isoflurane anesthesia in rats. Anesthesiology. 1992;77:920. [PubMed]
20. Van Gils M, Rosenfalck A, White S, Prior P, Gade J, Senhadji L, Thomsen C, Ghosh IR, Longford RM, Jensen K. Signal processing in prolonged EEG recordings during intensive care. IEEE Eng Med Biol Mag. 1997;16:56–63. [PubMed]
21. Lipping T, Jäntti V, Yli-Hankala A, Hartikainen K. Adaptive segmentation of burst-suppression pattern in isoflurane and enflurane anesthesia. Int J Clin Monit Comput. 1995;12:161–7. [PubMed]
22. Leistritz L, Jäger H, Schelenz C, Witte H, Putsche P, Specht M, Reinhart K. New approaches for the detection and analysis of electroencephalographic burst-suppression patterns in patients under sedation. J Clin Monit Comput. 1999;15:357–67. [PubMed]
23. Särkelä M, Mustola S, Seppänen T, Koskinen M, Lepola P, Suominen K, Juvonen T, Tolvanen-Laakso H, Jäntti V. Automatic analysis and monitoring of burst suppression in anesthesia. Int J Clin Monit Comput. 2002;17:125–34. [PubMed]
24. Rampil IJ, et al. I653 and isoflurane produce similar dose-related changes in the electroencephalogram of pigs. Anesthesiology. 1988;69:298. [PubMed]
25. Ching S, Purdon PL, Vijayan S, Kopell NJ, Brown EN. A neurophysiological–metabolic model for burst suppression. Proc Natl Acad Sci USA. 2012;109:3095–100. [PubMed]
26. Nakashima K, Todd MM, Warner DS. The relation between cerebral metabolic rate and ischemic depolarization: a comparison of the effects of hypothermia, pentobarbital, and isoflurane. Anesthesiology. 1995;82:1199. [PubMed]
27. Michenfelder JD. The interdependency of cerebral functional and metabolic effects following massive doses of thiopental in the dog. Anesthesiology. 1974;41:231. [PubMed]
28. Newberg LA, et al. The cerebral metabolic effects of isoflurane at and above concentrations that suppress cortical electrical activity. Anesthesiology. 1983;59:23. [PubMed]
29. Bruhn J, Bouillon TW, Shafer SL. Onset of propofol-induced burst suppression may be correctly detected as deepening of anaesthesia by approximate entropy but not by bispectral index. Br J Anaesth. 2001;87:505–7. [PubMed]
30. Zhang D, Jia X, Ding H, Ye D, Thakor NV. Application of tsallis entropy to EEG: quantifying the presence of burst suppression after asphyxial cardiac arrest in rats. IEEE Trans Biomed Eng. 2010;57:867–74. [PMC free article] [PubMed]
31. Löfhede J, Löfgren N, Thordstein M, Flisberg A, Kjellmer I, Lindecrantz K. Classification of burst and suppression in the neonatal electroencephalogram. J Neural Eng. 2008;5:402. [PubMed]
32. Smith AC, Brown EN. Estimating a state-space model from point process observations. Neural Comput. 2003;15:965–91. [PubMed]
33. Smith AC, Frank LM, Wirth S, Yanike M, Hu D, Kubota Y, Graybiel AM, Suzuki WA, Brown EN. Dynamic analysis of learning in behavioral experiments. J Neurosci. 2004;24:447–61. [PubMed]
34. Smith AC, Stefani MR, Moghaddam B, Brown EN. Analysis and design of behavioral experiments to characterize population learning. J Neurophysiol. 2005;93:1776–92. [PubMed]
35. Pawitan Y. In All Likelihood. Oxford: Oxford University Press; 2001.
36. Brown EN, Frank LM, Tang D, Quirk MC, Wilson MA. A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. J Neurosci. 1998;18:7411–25. [PubMed]
37. De Jong P, Mackinnon MJ. Covariances for smoothed estimates in state space models. Biometrika. 1988;75:601–2.
38. Meuret P, Backman SB, Bonhomme V, Plourde G, Fiset P. Physostigmine reverses propofol-induced unconsciousness and attenuation of the auditory steady state response and bispectral index in human volunteers. Anesthesiology. 2000;93:708. [PubMed]
39. Funk W, Hollnberger H, Geroldinger J. Physostigmine and anaesthesia emergence delirium in preschool children: a randomized blinded trial. Eur J Anaesthesiol. 2008;25:37–42. [PubMed]
40. Fan J, Gijbels I. Local Polynomial Modelling and its Applications. Vol. 66. London: Chapman and Hall; 1996.
41. Brown RG, Hwang PYC. Introduction to Random Signals and Applied Kalman Filtering: with MATLAB Exercises and Solutions. New York: Wiley; 1996.
42. Smith AC, et al. State-space algorithms for estimating spike rate functions. Comput Intell Neurosci. 2010;2010:426539. [PMC free article] [PubMed]
43. Liberman MY, Ching S, Chemail J, Brown EN. A closed-loop anesthesia delivery system for real-time control of burst suppression. J Neural Eng. 2013;10:046004. [PubMed]
44. Ching S, Westover BM, Liberman MY, Chemali JJ, Kenny J, Solt K, Purdon PL, Brown EN. Real-time closed loop control in a rodent model of medically-induced coma. Anesthesiology. 2013;87 at press. [PubMed]
45. Shanechi MM, Chemali JJ, Liberman MY, Solt K, Brown EN. A brain–machine interface for control of medically-induced coma. PLoS Comput Biol. 2013 at press. [PMC free article] [PubMed]