Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3770149

Formats

Article sections

- Abstract
- I. Introduction
- II. Model Definition
- III. Estimation Algorithm
- IV. Results
- V. Conclusion
- References

Authors

Related links

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

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

Jessica J. Chemali, Massachussets General Hospital ; Email: ude.tim.tatsoruen@40cjj;

The publisher's final edited version of this article is available at Conf Proc IEEE Eng Med Biol Soc

See other articles in PMC that cite the published article.

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 [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].

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.

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 *b _{i}* be the indicator of the
event at time

$$Pr({b}_{i}\mid {p}_{i},{x}_{i})={p}_{i}^{{b}_{i}}{(1-{p}_{i})}^{1-{b}_{i}},$$

(1)

where *p _{i}* is defined by the logistic
equation

$${p}_{i}=\frac{exp\left({x}_{i}\right)}{1+exp\left({x}_{i}\right)}.$$

(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 *p _{i}* is always between 0 and 1 as the
state

$${x}_{i}={x}_{i-1}+{\u220a}_{i},$$

(3)

where the *ε _{i}* are independent
Gaussian random variables with mean 0 and variance

Our objective is to estimate the state process *x _{i}*
and the parameter ${\sigma}_{e}^{2}$ for

$$\begin{array}{cc}\hfill p({b}_{1:N},x\mid {\sigma}_{e}^{2},{x}_{0})& ={\Pi}_{i=1}^{N}{p}_{i}^{{b}_{i}}{(1-{p}_{i})}^{1-{b}_{i}}{\left(\frac{1}{2\pi {\sigma}_{e}^{2}}\right)}^{\frac{1}{2}}\times exp(-\frac{1}{2}\frac{{({x}_{i}-{x}_{i-1})}^{2}}{{\sigma}_{e}^{2}})\hfill \\ \hfill & =p({b}_{1:N}\mid x)p(x\mid {\sigma}_{e}^{2},{x}_{0}).\hfill \end{array}$$

(4)

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
*b*_{1:N} and the estimates
*σ*^{2(l)} and ${x}_{0}^{\left(l\right)}$ of the parameters from iteration *l*:

$$Q({\sigma}_{e}^{2},{x}_{0}\mid {\sigma}^{2\left(l\right)},{x}_{0}^{\left(l\right)})=E[\text{log}[p({b}_{1:N},x\mid {b}_{1:N},{\sigma}^{2\left(l\right)},{x}_{0}^{\left(l\right)})=E({\Sigma}_{i=1}^{N}{b}_{i}{x}_{i}-\text{log}(1+\text{exp}\left({x}_{i}\right)\mid {b}_{1:N},{\sigma}^{2\left(l\right)},{x}_{0}^{\left(l\right)}))+E[{\Sigma}_{i=1}^{N}-\frac{1}{2}\frac{{({x}_{i}-{x}_{i-1})}^{2}}{{\sigma}_{e}^{2}}-\frac{N+1}{2}\text{log}\left(2\pi \right)-\frac{N+1}{2}\text{log}\left({\sigma}_{e}^{2}\right)-\frac{{x}_{0}^{2}}{2{\sigma}_{e}^{2}}\mid {b}_{1:N},{\sigma}_{e}^{2\left(l\right)},{x}_{0}^{\left(l\right)}].$$

(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*,

$${x}_{i\mid N}\equiv E[{x}_{i}\mid {b}_{1:N},{\sigma}_{e}^{2\left(l\right)},{x}_{0}^{\left(l\right)}],$$

(6)

and the covariances given the data up to time
*N*,

$${W}_{i\mid N}\equiv E[{x}_{i}^{2}\mid {b}_{1:N},{\sigma}_{e}^{2\left(l\right)},{x}_{0}^{\left(l\right)}]$$

(7)

and,

$${W}_{i,i-1\mid N}\equiv E[{x}_{i}{x}_{i-1}\mid {b}_{1:N},{\sigma}_{e}^{2\left(l\right)},{x}_{0}^{\left(l\right)}].$$

(8)

The expectation step consists of three parts. First we compute the
estimates of *x*_{i|i}
and ${\sigma}_{i\mid i}^{2}$ from the forward filter. Second, we use the backward filter to
compute *x*_{i|N} and ${\sigma}_{i\mid N}^{2}$. Finally, we use the state space covariance algorithm to
compute the covariances *W _{i|N}* and

Given the parameters estimates from iteration *l*,
this step estimates *x _{i|i}* and ${\sigma}_{i\mid i}^{2}$, meaning that it will estimate the state and the variance
at

$${x}_{i\mid i-1}={x}_{i-1\mid i-1},$$

(9)

$${\sigma}_{i\mid i-1}^{2}={\sigma}_{i-1\mid i-1}^{2}+{\sigma}_{e}^{2\left(l\right)}.$$

(10)

The posterior mode and variance are given by

$${x}_{i\mid i}={x}_{i\mid i-1}+{\sigma}_{i\mid i-1}^{2}({b}_{i}-{p}_{i\mid i}),$$

(11)

$${\sigma}_{i\mid i}^{2}={[{\left({\sigma}_{i\mid i-1}^{2}\right)}^{-1}+{p}_{i\mid i}(1-{p}_{i\mid i})]}^{-1}.$$

(12)

The initial conditions are ${x}_{0\mid 0}={x}_{0}^{\left(l\right)}$ and ${\sigma}_{0\mid 0}^{}={\sigma}^{2\left(l\right)}$. *p _{i|i}* corresponds to the mode
of the posterior distribution which we estimate recursively. This filter is
non-linear because

Using the posterior estimates, a fixed-interval smoothing algorithm
(FIS) [13][14] gives us the estimates
*x _{i|N}* for

$${x}_{i\mid N}={x}_{i\mid i}+{A}_{i}({x}_{i+1\mid N}-{x}_{i+1\mid i}),$$

(13)

$${A}_{i}={\sigma}_{i\mid i}^{2}{\left({\sigma}_{i+1\mid i}^{2}\right)}^{-1},$$

(14)

$${\sigma}_{i\mid N}^{2}={\sigma}_{i\mid i}^{2}+{A}_{i}^{2}({\sigma}_{i+1\mid N}^{2}-{\sigma}_{i+1\mid i}^{2}).$$

(15)

The initial conditions are *x _{N|N}* and ${\sigma}_{N\mid N}^{2}$ previously estimated in the filter algorithm.

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

$${\sigma}_{i,j\mid N}={A}_{i}{\sigma}_{i+1,j\mid N},$$

(16)

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

$${W}_{i,i-1\mid N}={\sigma}_{i,i-1\mid N}+{x}_{i\mid N}{x}_{i-1\mid N}$$

(17)

and,

$${W}_{k,k\mid K}={\sigma}_{i\mid N}^{2}+{x}_{i\mid N}^{2}.$$

(18)

The maximization step is

In the maximization step, the complete data log likelihood is maximized
with respect to *σ*^{2} and
*x*_{0}. This yields the estimates

$${\sigma}^{2(l+1)}={(N+1)}^{-1}(2({\Sigma}_{i=2}^{N}{W}_{i\mid N}-{\Sigma}_{i=2}^{N}{W}_{i-1,i\mid N})+\frac{3}{2}{W}_{1\mid N}-{W}_{N\mid N})$$

(19)

and,

$${x}_{0}^{(l+1)}=\frac{1}{2}{x}_{1\mid N}.$$

(20)

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 *x*_{0} and ${\sigma}_{2}^{2}$ 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:

$$f(p\mid {x}_{i\mid N},{\sigma}_{i\mid N}^{2},{\sigma}_{i\mid N}^{2})={\left[{\left(2\pi {\sigma}_{i\mid N}^{2}\right)}^{0.5}p(1-p)\right]}^{-1}\times exp(-\frac{1}{2{\sigma}_{i\mid N}^{2}}{[log\left[p{(1-p)}^{-1}\right]-{x}_{i\mid N}]}^{2}).$$

(21)

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).

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*(*x _{i}* >

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*(*x _{i}*
<

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: ude.tim.tatsoruen@40cjj.

K.F. Kevin Wong, Massachussets General Hospital ; Email: ude.tim.tatsoruen@gnow.

Ken Solt, Harvard Medical School/Massachussets General Hospital ; Email: gro.srentrap@tlosk.

Emery N. Brown, Harvard Medical School/Massachussets General Hospital and professor of computational neuroscience and health sciences and technology at MIT ; Email: ude.tim.tatsoruen@bne.

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |