PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2013 September 11.
Published in final edited form as:
PMCID: PMC3770149
NIHMSID: NIHMS504743

A state-space model of the burst suppression ratio

Jessica J. Chemali, research technologist, K.F. Kevin Wong, post-doctoral fellow, Ken Solt, assistant professor of anesthesia, and Emery N. Brown, professor of anesthesia

Abstract

Burst suppression is an electroencephalogram pattern observed in states of severely reduced brain activity, such as general anesthesia, hypothermia and anoxic brain injuries. The burst suppression ratio (BSR), defined as the fraction of EEG spent in suppression per epoch, is the standard quantitative measure used to characterize burst suppression. We present a state space model to compute a dynamic estimate of the BSR as the instantaneous probability of suppression. We estimate the model using an approximate EM algorithm and illustrate its application in the analysis of rodent burst suppression recordings under general anesthesia. Our approach removes the need to artificially average the ratio over long epochs and allows us to make formal statistical comparisons of burst activity at different time points. Our state-space model suggests a more principled way to analyze this key EEG feature that may offer more informative assessments of its associated brain state.

I. Introduction

Burst suppression is an electroencephalogram (EEG) pattern consisting of alternating high amplitude bursts and low amplitude or isoelectric periods (see Fig.1). Burst suppression can be observed in deep state of general anesthesia, coma associated with diffuse brain injuries, induced hypothermia, epilepsy due to Ohtahara's syndrome [1], and postasphyctic newborn babies [2]. It is usually caused by physical damage (e.g brain injury, lesions), or reversible general anesthetic effects. In burst suppression less than 5% of the cortical neurons and less than 40% of the thalamic neurons are active during suppression [3].

Fig. 1
The EEG signals of a rat, awake, under 2.5% and 4 % Isoflurane.

Burst suppression patterns (BSPs) differ in the durations of the alternating periods. For instance, in the case of subacute sclerosing panencephalitis, the duration between bursts tends to shorten with the progression of the disease [4]. In the case of general anesthesia, BSPs produced by different anesthetics show clear differences [5]. The burst suppression ratio (BSR), a measure of the fraction of time spent in suppression per epoch is traditionally used to characterize BSPs. The conventions for computing BSR are chosen empirically. The voltage threshold, which discriminates between low and high amplitude can be set at 15 [6]-[8] or at 5 [9]-[11] μV , and the minimal suppression time that the EEG has to be isoelecric is set at 100 milliseconds in some cases [6]-[8] and 0.5 seconds in others [9]-[11]. Likewise, the length of the epoch over which the fraction of suppression is computed can take values from 4 [9]-[11] to 15 [6]-[8] seconds. Finally the number of epochs that is averaged is typically chosen to span 60 seconds [6]-[11] resulting in only one data point per minute. The ease of computation and its real-time application make this method appealing. However, it offers a coarse characterization of the EEG pattern and limits the information extracted from the data sets. For instance, when analyzing the modulation of the BSP caused by certain drugs and inferring from it the anesthetic state of the brain, important dynamics taking place over several seconds can be hidden by averaging. This approach artificially segments the EEG signal into epochs and averages sequences whereas the signal is continuous and changes dynamically. Also, because different conventions are used to estimate the BSR, comparing results across studies is challenging [12]. The need for a more precise definition of BSR is thus motivated by the interest in inferring the dynamics of brain states from the relevant variations in BSPs. A more rigorous definition also facilitates achieving consistency between results across different studies, opening the door for insightful comparisons.

The BSP can be described by a series of 1s (suppression) and 0s (no suppression). This definition casts the problem in a dynamic Bernoulli process framework. Therefore, we can analyze it using a state-space model to characterize the burst suppression ratio as the probability of suppression as a function of time. A Bernoulli probability model is used to describe the suppression occurrences and a Gaussian state equation describes the unobservable brain state whose evolution we want to track. We estimate the model using an approximate expectation maximization (EM) algorithm [13][14], and compute the dynamic estimate of the BSR. We illustrate our method by computing the BSR before and after the administration of the cholinergic agonist physostigmine. We show how our approach may be used to conduct formal statistical analyses of the BSR.

II. Model Definition

A state space model is completely characterized by two equations: a state equation and an observation equation. The state equation defines the unobservable state process whose evolution we wish to track over time. The state in our case represents the brain state of the subject. We define our state to be positively correlated with the BSR. That is as the state increases the BSR increases and as the state decreases the BSR decreases. The observation equation describes how the observed data are related to the unobservable state. In our case, the observed data are the sequence of 1s and 0s which represent the suppression or not of a burst respectively. We define a burst suppression to be a sequence of low amplitude EEG voltages over at least a 100 milliseconds. Our objective is thus to estimate the states and the parameters of the model from the observed data and use these estimates to compute the dynamic estimate of the BSR with confidence intervals. Our experiment consists of N discrete time steps indexed 1, ..., N. Let bi be the indicator of the event at time i which is 1 if the event is suppressed and 0 otherwise. Let pi denote the probability of a suppression at time i. In our state space model we assume that the probability pi is governed by an unobservable state variable xi. The observation model can be described by the Bernoulli probability mass function

equation M1
(1)

where pi is defined by the logistic equation

equation M2
(2)

The logistic function is commonly used to model how the probability of an event is affected by one or many explanatory variables. In our case it ensures that the probability pi is always between 0 and 1 as the state xi varies over all real numbers. We define the unobserved state model by a random walk

equation M3
(3)

where the εi are independent Gaussian random variables with mean 0 and variance σ2. This definition of the state model provides a stochastic continuity constraint, which insures that states that are close in time are close in value. In our case this guarantees continuity of the burst probabilities. The parameter σ2 governs how rapidly changes can occur between adjacent states. To estimate the set of states and σ2 we use an expectation maximization (EM) algorithm for point processes and binary time series [13][14].

III. Estimation Algorithm

Our objective is to estimate the state process xi and the parameter equation M4 for i = 1, ..., N. This will allow us to compute the BSR estimate with its confidence intervals. We treat the initial condition x0 as a parameter. The states and parameter estimates are computed by maximizing the expectation of the complete data log likelihood. The complete data likelihood is

equation M5
(4)

The EM algorithm iterates between an expectation or E-step and a maximization or M-step. The expectation step is

A. Expectation Step

At iteration l + 1 we compute in the expectation step the expected value of the complete data log likelihood given the complete data b1:N and the estimates σ2(l) and equation M6 of the parameters from iteration l:

equation M7
(5)

We can see upon further expanding the right hand side that we need to estimate three quantities for i = 1, ..., N

The expectation of the state variable given the data up to time N,

equation M8
(6)

and the covariances given the data up to time N,

equation M9
(7)

and,

equation M10
(8)

The expectation step consists of three parts. First we compute the estimates of xi|i and equation M11 from the forward filter. Second, we use the backward filter to compute xi|N and equation M12. Finally, we use the state space covariance algorithm to compute the covariances Wi|N and Wi,i–1|N.

1) Forward Filter

Given the parameters estimates from iteration l, this step estimates xi|i and equation M13, meaning that it will estimate the state and the variance at i looking at data from the start of the experiment up to i using a non linear forward filter algorithm [13][14]. The one step prediction mean and variance are given by

equation M14
(9)
equation M15
(10)

The posterior mode and variance are given by

equation M16
(11)
equation M17
(12)

The initial conditions are equation M18 and equation M19. pi|i corresponds to the mode of the posterior distribution which we estimate recursively. This filter is non-linear because xi|i appears on both sides of (11).

2) Backward Filter

Using the posterior estimates, a fixed-interval smoothing algorithm (FIS) [13][14] gives us the estimates xi|N for i = N – 1, ..., 1 , meaning that the estimate at time i is based on all the data up to time N. Our final estimate of the state will be thus a Gaussian distributed variable with mean xi|N and variance equation M20. The FIS is

equation M21
(13)
equation M22
(14)
equation M23
(15)

The initial conditions are xN|N and equation M24 previously estimated in the filter algorithm.

3) State Space Covariance Algorithm

Finally, the state space covariance algorithm [15] gives σi,j|N:

equation M25
(16)

where 1 < i < j < N. The covariances are thus given by

equation M26
(17)

and,

equation M27
(18)

The maximization step is

B. Maximization Step

In the maximization step, the complete data log likelihood is maximized with respect to σ2 and x0. This yields the estimates

equation M28
(19)

and,

equation M29
(20)

Dynamic BSR Estimate

The algorithm iterates between the Expectation and Maximization steps until convergence. The convergence criteria we use are the same as those developed by [12]. The fixed-interval smoothing algorithm evaluated at the maximum likelihood estimates x0 and equation M30 together with the logistic equation (2) give us the probability of a suppression at time i for i = 1, ..., N. From there we can compute the probability density function which corresponds to the dynamic BSR estimate and is given by:

equation M31
(21)

IV. Results

Fig. 2 illustrates the results of our algorithm. The data used were recorded from a rat under stable 2% isoflurane for sixty minutes. A four-minute window centered at the time of physostigmine injection is shown in Fig. 2. We adapted the algorithm of Vijn and colleagues [6] to extract the burst suppression binary sequence by setting the minimum burst suppression window and the epoch size to 100 milliseconds. This binary sequence is the input to our algorithm. The analysis was conducted using the Matlab software (MathWorks, Natick, MA) available at our website (http://neurostat.mgh.harvard.edu/BehavioralLearning/Matlabcode).

Fig. 2
Panel A shows the estimated BSR (red curve) with its 95% confidence intervals (black curves). The BSR declines appreciably after administration of physostigmine at time 0. Panel B shows a scatter plot of the significant differences (see text) between ...

Panel A in Fig. 2 shows the dynamic BSR curve with its 95% confidence intervals. Prior to the administration of physostigmine at time 0, the BSR oscillates approximately between 0.2 and 0.8 with a median value of approximately 0.6. This suggests that the EEG spent the majority of its time in the suppressed state. After time 0, there is a steep decrease in the BSR. From this point on the BSR oscillates between 0.05 and 0.10. This is consistent with a decrease in suppression or equivalently, an increase in arousal. Because we estimated the joint distribution of the state process we can evaluate Pr(xi > xj) for any 0 ≤ i < jN which is equivalent to the probability that the BSR at time i is greater than the BSR at time j because the transformation between the state variable xi and the probability of a suppression pi is monotonic. Panel B in Fig.2 is a scatter plot showing the pairwise significant differences between the BSR across time. The red and black dots correspond to Pr(xi < xj) = 0.975 and Pr(xi > xj) = 0.975 respectively, where i is the index across the xaxis and j is the index across the yaxis.

To illustrate how to read Fig. 2B, notice that the BSR decreases at time −0.75 min. The red vertical strip at time −0.75 min show that the BSR at all times from −2 min to −0.75 min are significantly greater than the BSR at −0.75 min to −0.80 min. In contrast, the black vertical strip at the −0.5 min shows that all BSR prior to this time are significantly smaller. This observation is consistent with the BSR estimates in Fig. 2A. The red dots clearly show that Pr(xi < xj) = 0.975 for i ≥ 0 and j < 0, which means that the BSR after the injection of the physostigmine is significantly smaller than before the injection. This abrupt drop in the BSR is expected because physostigmine is a cholinergic agonist that induces increased arousal.

V. Conclusion

We have presented a state space model in which we define BSR dynamically as the probability of a suppression period as a function of time. We estimated the model using an approximate EM algorithm and illustrated its application in the analysis of rodent burst suppression recordings under general anesthesia. Our approach removes the need to average the ratio over long epochs. It also provides a continuous estimate of the probability of a burst and the joint distribution of the burst suppression states. The latter allows us to make formal statistical comparisons between burst activity at any time points. Our state-space model suggests a more principled way to analyze this key EEG feature that may offer more informative assessment of this important brain state.

Acknowledgments

This work was supported by the Massachusetts General Hospital Department of Anesthesia, Critical Care and Pain Medicine and by NIH grants DP1 OD003646-01 and R01 MH071847 to E.N.B., and NIH grant K08 GM094394 to K.S.

Contributor Information

Jessica J. Chemali, Massachussets General Hospital ; jjc04/at/neurostat.mit.edu.

K.F. Kevin Wong, Massachussets General Hospital ; wong/at/neurostat.mit.edu.

Ken Solt, Harvard Medical School/Massachussets General Hospital ; ksolt/at/partners.org.

Emery N. Brown, Harvard Medical School/Massachussets General Hospital and professor of computational neuroscience and health sciences and technology at MIT ; enb/at/neurostat.mit.edu.

References

1. Brown EN, Lydic R, Schiff ND. General anesthesia, sleep, and coma. N. Engl. J. Med. 2010;363(27):2638–2650. [PMC free article] [PubMed]
2. Lofhede J. PhD. dissertation. Dept. Signals and Systems, Chalmers University of Technology; Goteborg, Sweden: 2009. The EEG of the neonatal brain classification of background activity.
3. Steriade M, Amzica F, Contreras D. Cortical and thalamic cellular correlates of EEG burst suppression,” Electroencephalogr. Clin. Neurophysiol. 1994;90(1):1–16. [PubMed]
4. Boggs JG. Generalized EEG Waveforms Abnormalities. http://emedicine.medscape.com/article/1140075-overview, Mar. 12, 2009 [Apr.14, 2011]
5. Akrawi WP, William JC, Drummond CJ, Kalkman, Patel PM. A Comparison of the Electrophysiologic Characteristics of EEG Burst-Suppression as Produced by Isoflurane, Thiopental, Etomidate, and Propofol. J. Neurosurg. Anesthesiol. 1996;8:40–46. [PubMed]
6. Vijn PC, Sneyd JR. I.v. anaesthesia and EEG burst suppression in rats: bolus injections and closed-loop infusions. Br. J. Anaesth. 1998;81:415–421. [PubMed]
7. Luo T, Leung LS. Basal forebrain histaminergic transmission modulates electroencephalographic activity and emergence from isoflurane anesthesia. Anesthesiology. 2009;111:725–733. [PubMed]
8. Van Den Broek PL, Van Rijn CM, Van Egmond J, Coenen AM, Booij LH. An effective correlation dimension and burst suppression ratio of the EEG in rat. correlation with sevoflurane induced anaesthetic depth. Eur. J. Anaesthesiol. 2006;23:391–402. [PubMed]
9. Rampil IJ. A primer for EEG signal processing in anaesthesia. Anesthesiology. 1998;89:980–1002. [PubMed]
10. Rampil IJ, Laster MJ. No correlation between quantitative electroencephalographic measurements and movement response to noxious stimuli during isoflurane anesthesia in rats. Anesthesiology. 77:920–925. [PubMed]
11. Doyle PW, Matta BF. Burst suppression or isoelectric encephalogram for cerebral protection: evidence from metabolic suppression studies. Br. J. Anaesth. 1999;83(4):580–584. [PubMed]
12. Van De Velde M. PhD. dissertation. Eindhoven University of Technology; Eindhoven, Netherlands: 2000. Signal validation in electroencephalography research.
13. Smith AC, Brown EN. Estimating a state-space model from point process observations. Neural. Comp. 2003;15:965–991. [PubMed]
14. 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–461. [PubMed]
15. De Jong P, Mackinnon MJ. Covariances for smoothed estimates in state space models. Biometrika. 1988;75:601–602.