|Home | About | Journals | Submit | Contact Us | Français|
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.
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 , and postasphyctic newborn babies . 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 .
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 . In the case of general anesthesia, BSPs produced by different anesthetics show clear differences . 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 - or at 5 - μV , and the minimal suppression time that the EEG has to be isoelecric is set at 100 milliseconds in some cases - and 0.5 seconds in others -. Likewise, the length of the epoch over which the fraction of suppression is computed can take values from 4 - to 15 - seconds. Finally the number of epochs that is averaged is typically chosen to span 60 seconds - 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 . 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 , 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.
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
where pi is defined by the logistic equation
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
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 .
Our objective is to estimate the state process xi and the parameter 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
The EM algorithm iterates between an expectation or E-step and a maximization or M-step. The expectation step is
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 of the parameters from iteration l:
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,
and the covariances given the data up to time N,
The expectation step consists of three parts. First we compute the estimates of xi|i and from the forward filter. Second, we use the backward filter to compute xi|N and . Finally, we use the state space covariance algorithm to compute the covariances Wi|N and Wi,i–1|N.
Given the parameters estimates from iteration l, this step estimates xi|i and , 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 . 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 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).
Using the posterior estimates, a fixed-interval smoothing algorithm (FIS)  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 . The FIS is
The initial conditions are xN|N and previously estimated in the filter algorithm.
Finally, the state space covariance algorithm  gives σi,j|N:
where 1 < i < j < N. The covariances are thus given by
The maximization step is
In the maximization step, the complete data log likelihood is maximized with respect to σ2 and x0. This yields the estimates
The algorithm iterates between the Expectation and Maximization steps until convergence. The convergence criteria we use are the same as those developed by . The fixed-interval smoothing algorithm evaluated at the maximum likelihood estimates x0 and 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:
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  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).
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 < j ≤ N 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 x – axis and j is the index across the y – axis.
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.
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.
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.
Jessica J. Chemali, Massachussets General Hospital ; Email: jjc04/at/neurostat.mit.edu.
K.F. Kevin Wong, Massachussets General Hospital ; Email: wong/at/neurostat.mit.edu.
Ken Solt, Harvard Medical School/Massachussets General Hospital ; Email: 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 ; Email: enb/at/neurostat.mit.edu.