|Home | About | Journals | Submit | Contact Us | Français|
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.
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.
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.
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 . Observed in profound general anesthesia , coma , hypothermia , epilepsy due to Ohtahara’s syndrome , and postasphyctic newborns , 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 . Figures 1(A) and (C) show respectively a 5-min and a 1-min segment of burst suppression induced by the anesthetic propofol.
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 . 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 . 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 . Persistence of burst suppression patterns in the EEG of patients in coma is commonly associated with a poor prognosis . 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 . 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 [18–23]. 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 . 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) . 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, 26–28]. Other approaches to characterizing burst suppression have included entropy measures [29, 30] and artificial neural networks and support vector machines .
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 [32–34]. 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.
To formulate our state-space model, we follow the state-space paradigm for analyzing point processes and binary time-series in [32–34]. 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 Δ = nδ 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
where pi, the BSP, is
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
where the εi are independent, Gaussian random variables with mean 0 and variance . 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 governs how rapidly the BSP can change; the larger (smaller) the value of the more rapidly (slowly) the state and the BSP can change.
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
given the initial state x0. The joint probability mass function of the observed suppression events given the states is
Our objective is to estimate using maximum likelihood (ML) the state process x and the parameters 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 [32–34]. 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 . The complete data likelihood is
The EM algorithm is
At iteration + 1 we compute in the expectation step the expected value of the complete data log likelihood given the data b and the estimates and of the parameters from iteration :
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
and the covariances of the state variables conditioned on the data up to time I are
To compute these quantities efficiently, we divide the expectation step into three parts. First we compute the estimates of xi|i and using a binary filter algorithm. Second, we use a fixed interval smoothing (FIS) algorithm to compute xi|I and . Finally, we use the state-space covariance algorithm to compute the covariances Wi|I and Wi,i−1|I.
The one step prediction mean and variance are given by
The posterior mode and variance are given by
The initial conditions are 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.
Using the estimates from the binary filter, a FIS algorithm [32, 33] gives the state and variance estimates xi|I and 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 . The FIS is
The initial conditions are xI|I and , computed at the final step of the binary filter algorithm.
Finally, using the state-space covariance algorithm  we compute σi, j|I :
To carry out the maximization step at iteration + 1 we let and assume that τ has a gamma prior density defined as
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
We chose the gamma distribution for the prior because it is flexible and because it is the conjugate distribution for τ in the Gaussian likelihood . 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
We maximize (22) with respect to x0 and find
The algorithm iterates between the expectation and maximization steps until convergence. The convergence criteria we use are the same as those developed by . The ML estimates of τ, or equivalently , and x0 are respectively and , 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 evaluated at their ML estimates, we term them respectively the BSP filter and smoothing algorithms.
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
The 95% confidence intervals are computed from the cumulative density of (25) by identifying its 2.5th and 97.5th percentiles.
Given the ML estimates of x0 and , 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 . 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 . 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
We can then draw M samples from the Gaussian distribution with mean
and covariance matrix
and count the number G of instances in which xj|I > xi|I. The estimate of the probability of interest is
In our analysis we chose M=10 000.
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.
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.
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 . 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.
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
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.
Physostigmine is a cholimimetic drug that has been used in anesthesiology research to induce emergence from general anesthesia  and in anesthesiology practice to treat emergence delirium , 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 . 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).
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.
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 , 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 < j ≤ I. 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.
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 . 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)).
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  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 , 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 .
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 governs the degree of smoothing in the BSP estimates. Larger (smaller) values of allow for less (more) smoothness in the BSP time course. Placing a prior distribution on 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 , 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, bi − npi|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 , which is the gain in the BSP filter, governs how much the innovation is weighted in computing the new update. Because 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 . These local adaptive features of the BSP algorithms, which are characteristic of Kalman filter-like algorithms , 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 , and as well as or better than more elaborate smoothing algorithms that have an automatic bandwidth selection criterion .
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 . 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  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 . 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 [43–45] will be the topics of future reports.
This work was supported by: NIH Awards DP1-OD003646 (ENB), DP2-OD006454 (PP), K08-GM094394 (KS), and the Burroughs-Wellcome Fund Award 1010625 (SC).
Some figures may appear in colour only in the online journal