Emergence of spatially heterogeneous burst suppression in a neural field model of electrocortical activity.

Burst suppression in the electroencephalogram (EEG) is a well-described phenomenon that occurs during deep anesthesia, as well as in a variety of congenital and acquired brain insults. Classically it is thought of as spatially synchronous, quasi-periodic bursts of high amplitude EEG separated by low amplitude activity. However, its characterization as a "global brain state" has been challenged by recent results obtained with intracranial electrocortigraphy. Not only does it appear that burst suppression activity is highly asynchronous across cortex, but also that it may occur in isolated regions of circumscribed spatial extent. Here we outline a realistic neural field model for burst suppression by adding a slow process of synaptic resource depletion and recovery, which is able to reproduce qualitatively the empirically observed features during general anesthesia at the whole cortex level. Simulations reveal heterogeneous bursting over the model cortex and complex spatiotemporal dynamics during simulated anesthetic action, and provide forward predictions of neuroimaging signals for subsequent empirical comparisons and more detailed characterization. Because burst suppression corresponds to a dynamical end-point of brain activity, theoretically accounting for its spatiotemporal emergence will vitally contribute to efforts aimed at clarifying whether a common physiological trajectory is induced by the actions of general anesthetic agents. We have taken a first step in this direction by showing that a neural field model can qualitatively match recent experimental data that indicate spatial differentiation of burst suppression activity across cortex.

1. Introduction

Over the many years since its discovery in humans (Berger, 1929, 1930; Adrian and Matthews, 1934), the electroencephalogram (EEG) has been shown to be a sensitive, and often specific, indicator of brain state and function (Schomer and Lopes da Silva, 2010). In the case of the deeply inactivated brain, whether through trauma or medical intervention, a burst suppression pattern is typically observed (Niedermeyer, 2009; Ching et al., 2012). Consisting of quasi-periodic alternations of high amplitude periods of spiking activity with low amplitude periods that are near isoelectric, the burst suppression pattern is associated with a range of central insults or interventions that include cortical deafferentation (Henry and Scoville, 1952; Kellaway et al., 1966; Lukatch and MacIver, 1996), cerebral ischaemia (Bauer et al., 2013), deep coma (Young, 2000), various infantile encephalopathies (Grigg-Damberger et al., 1989), the final stages of deteriorated status epilepticus (Treiman et al., 1990), hypothermia (Stecker et al., 2001), and high levels of many anesthetic and sedative drugs (Schwartz et al., 1989; Akrawi et al., 1996).

The burst suppression pattern can show a significant degree of variation depending on its aetiology. For example, in the case of infantile hypoxic-ischemic encephalopathy the burst suppression pattern can be quite complex; and due to significant variability in the amplitude of individual bursts a clear transition to suppression may not readily be apparent (Lamblin et al., 2013). In contrast, in deep anesthesia bursts are typically separated by clear isoelectric periods, the duration (relative and absolute) of which increases systematically with increasing anesthetic level. This systematic dependence on anesthetic level can be utilized in the treatment of status epilepticus (Kalviainen et al., 2005) and the management of brain trauma in the intensive care setting (Doyle and Matta, 1999) by defining an endpoint in which more than 50% of an EEG recording consists of suppressions.

In what follows we will provide an overview of the phenomenon of burst suppression and summarize the current understanding regarding its physiological genesis. This will then be followed by a detailed outline of a neural field model developed to describe the emergence of burst suppression during anesthesia, which notably, and for the first time, incorporates the empirically realistic modeling of a general anesthetic agent (isoflurane) and the spatio-temporal propagation of cortical activity.

1.1. Physiological basis of burst suppression

Despite its clear aetiological associations and clinical utility, little is known about the physiological mechanisms responsible for the genesis of burst suppression (Liley and Walsh, 2013). On the basis of brain slice and in vivo animal studies, a number of hypotheses have been advanced with sometimes contradictory conclusions. For example both increases (Steriade et al., 1994) and decreases (Ferron et al., 2009) in GABAergic inhibitory activity have been speculated to have causal roles in the onset of burst suppression. Supporting reductions in inhibition are in vivo whole-brain animal studies suggesting that enhanced network excitability (Detsch et al., 2002; Hudetz and Imas, 2007; Kroeger and Amzica, 2007; Land et al., 2012), possibly mediated through alterations in extracellular calcium (Kroeger and Amzica, 2007), is responsible for driving transitions between low amplitude quiescence and high amplitude bursting. The study of Land et al. (2012) is particularly relevant in this regard. Not only do they report that auditory and visual stimuli readily evoke burst activity in visual cortex (V1) and subiculum during deep anesthesia in rats, but (i) such excitability does not occur in the absence of burst suppression, (ii) V1 and subiculum bursting, in response to the cortically remote auditory stimulus, emerges abruptly with increasing anesthetic (isoflurane) concentration, and (iii) hysteresis occurs in both stimulus-induced and spontaneous bursting during isoflurane wash-in and wash-out. Thus, the phenomenon of burst suppression might be explicable in terms of the emergence of propagating excitability through a dynamical bifurcation parametrically regulated by isoflurane concentration.

Clinically it is well-established that bursting responses during burst suppression in deep anesthesia can be readily evoked by noxious and sensory stimulation, thus further implicating a role for alterations in cortical excitability in the genesis of burst suppression. For example, in probably the first systematic study on evoked bursts, Yli-Hankala et al. (1993) observed that a vibratory stimulus applied to the palm of the hand was readily able to evoke electroencephalographic bursts in patients during moderately deep isoflurane anesthesia. Subsequently it has been found that a range of visual, auditory, tactile and noxious stimuli are able to evoke electroencephalographic bursts during deep anesthesia in which burst suppression has been variously induced with isoflurane (Hartikainen et al., 1995), sevoflurane (Jantti et al., 1998) or propofol (Huotari et al., 2004).

Complementing this empirical and clinical research are recent modeling studies (Ching et al., 2012; Liley and Walsh, 2013), which suggest that the onset or unmasking of slow and activity-dependent modulations of network excitability might account for burst suppression patterns. Because it is observed that the spectral characteristics of the EEG just prior to the onset of the anesthesia-induced burst suppression are conserved in the bursts, such theoretical approaches typically modulate the oscillatory system that accounts for the dynamical emergence of the resting and anesthetic EEG. In order to simulate burst suppression during deep propofol anesthesia, Ching et al. (2012) utilize a thalamo-cortical model based on individual neurons previously developed to account for the propofol-induced emergence of frontal alpha-spindle activity (Ching et al., 2010). This model is then augmented with a slow adenosine triphosphate (ATP) gated potassium membrane current, which is hence regulated by the activity-dependent metabolic production rate of ATP. By assuming that propofol down-regulates neuronal firing through enhanced synaptic inhibition, thus leading to an autoregulatory decrease in cerebral metabolism and hence ATP production, the modulatory effect of this potassium current is magnified such that bursting emerges. In contradistinction to this model, in which bursting arises due to essentially intrinsic changes of neuronal excitability, Liley and Walsh (2013) developed a model in which bursting arises as a consequence of the slow activity-dependent modulation of synaptic efficacy. In this model the effects of synaptic resource depletion (receptor desensitization and synaptic vesicle depletion) and recovery during periods of sustained neuronal population activity act to slowly modulate neuronal population excitability. This mechanism comes to the fore in anesthesia because the general reduction of cortical activity allows the synaptic neurotransmitter reservoirs to fill up, potentiating excitation until it is sufficiently strong to induce feedback bursts of excitation, followed by suppression as the thereby depleted reservoirs refill. In support of this mechanism are the activity-dependent alterations in synaptic efficacy that have been been observed in vivo in recordings in cats during burst suppression induced with isoflurane (Kroeger and Amzica, 2007).

1.2. Spatio-temporal features of burst suppression

Because burst suppression is classically characterized as being a spatially homogeneous phenomenon (Brenner, 1985; An et al., 1996; Lewis et al., 2013), on the basis of near simultaneous burst onset and offset across scalp electrode derivations, little attention has been paid to its spatio-temporal features until recently. Motivated by the inability of scalp electroencephalography to reveal the fine structure of cortical dynamics, due to the spatial blurring induced by volume conduction, Lewis et al. (2013) chose to investigate the spatiotemporal features of burst suppression using intracranial electrocortigraphy (ECoG) in medically intractable epilepsy patients. Five patients, implanted with a range of subdural strip, grid and depth electrodes as part of a standard clinical monitoring procedure, had recordings collected throughout the induction of anesthesia with propofol, during explantation surgery. Burst onset and offset was observed to be visibly asynchronous across recording electrodes, with the absolute difference in burst onset time in general an increasing function of inter-electrode distance. Interestingly, not all recording electrodes would participate in such asynchronous bursting. It was found that burst onsets were visibly clustered across channels such that bursting could either be confined to a small subset of nearby electrodes (“local” bursting) or spread to involve the whole electrode grid (“global” bursting), with more distantly separated electrode pairs less likely to share a burst (based on burst onset within some time window).

It has been speculated that the appearance of spatially inhomogeneous bursting might be a reflection of the differential sensitivity of specific thalamo-cortical networks to anesthetic action. However, another possibility is that the spatially heterogeneous nature of this bursting arises as a feature of the axonal propagation of activity through cortex. In support of such a speculation is the developmental emergence of isoflurane-induced burst suppression in rats. It is conjectured that it is the development of short- and long-range horizontal connections between pyramidal neurons in infra-granular cortical layers, which is the critical factor in determining the appearance of isoflurane-induced burst suppression in the second postnatal week (Sitdikova et al., 2014). Further implicating the role that altered propagation may have in determining the physiological features of anesthetic action are reports that document the effects anesthetics have on nerve conduction—both centrally and peripherally. While peripherally it is generally assumed that anesthetics principally depress spinal motoneuron excitability, as assessed by reductions in F-wave amplitudes (Friedman et al., 1996; Rampil and King, 1996), there are a number of reports documenting the significant effects of anesthetic agents in either increasing (cyclopropane, nitrous oxide, diethyl ether) (Rosner et al., 1971) or reducing (pentobarbital, desflurane, enflurane, halothane) (Rampil and King, 1996; Oh et al., 2010; Nowicki et al., 2013) nerve conduction velocity at clinical levels, as assessed by increases in F-wave latency. Centrally, there is some evidence that volatile anesthetics may preferentially depress nerve conduction in unmyelinated axons (Berg-Johnsen and Langmoen, 1986; Mikulec et al., 1998). For instance, isoflurane was found to induce a conduction block in 20–30% of the unmyelinated fibers in the CA1 region of the rat hippocampus at clinical concentrations, as well as having a 1% effect on the actual conduction velocity (Berg-Johnsen and Langmoen, 1986). On the basis of empirical evidence indicating that the cortico-cortical fiber system is comprised of a mixture of myelinated and unmyelinated fibers, cf. Bojak and Liley (2010) and references therein, we hence expect mean cortical axonal conduction velocity to increase slightly, due to the reduction in the proportion of low conduction velocity unmyelinated fibers, but nevertheless anticipate cortico-cortical synaptic connectivity to be attenuated.

1.3. Necessity of large-scale cortical models

Regardless of the specific changes in cortical axonal conduction induced by anesthetics, it is clear that any theoretical attempt to account for burst suppression and its spatial inhomogeneity must explicitly incorporate the spatial extent of cortex. While the construction of a biophysically-based neuronal network model might seem an obvious starting point, numerical tractability and parametric uncertainties militates against the utility of this approach both from a descriptive and an explanatory perspective. For example, to meaningfully accommodate the extent of the spatially heterogeneous burst suppression seen in Lewis et al. (2013), we would need to model ~ 109 neurons and ~ 1012 synapses. While computations at this scale may be at the edge of feasibility for the largest supercomputers, we cannot reasonably expect such massive computations to be used for all the myriad specific research agendas in computational neuroscience any time soon. Even if such resources were readily available we would still be unable to specify the microcircuitry realistically at this level of detail for such a sizeable part of cortex. Such a problem will persist even if our computational capabilities continue to grow exponentially.

Fortunately, by considering the behavior of populations of neurons at mesoscopic scales, a variety of numerically tractable modeling approaches can be motivated physiologically and anatomically, cf. the reviews of Deco et al. (2008), Coombes (2010), Liley et al. (2012), and Liley (2013). These neural population models, referred to as neural mass models if localized and neural field or mean field models if spatially continuous and extensive, usually aim to describe the dynamical evolution of mean quantities (such as soma membrane potential or firing rate) defined over some suitable spatial domain or scale. Because these models average the activity of many thousands of neurons, they are well-suited as frameworks for understanding the meso- and macroscopic neural activity recorded, or inferred by, ECoG, EEG, magnetoencephalography (MEG) and the blood-oxygen level dependent (BOLD) contrast of functional magnetic resonance imaging (fMRI) (Bojak and Breakspear, 2013). Since the pioneering work of Walter Freeman (Freeman, 1975), this approach has flourished and has resulted in a number of important neural field models aimed at explaining the dynamical genesis of the mammalian EEG (Wilson and Cowan, 1973; Lopes da Silva et al., 1974; Nunez, 1974; Liley et al., 2002; Robinson et al., 2004). Broadly speaking, all these models are able to generate oscillatory activity through reverberant feedforward and feedback synaptic activity between excitatory and inhibitory neuronal populations. We choose to utilize the neural field model of Liley et al. (2002) as a framework for better understanding the spatial heterogeneity of bursting during anesthesia because (i) it has been previously employed to account for a number of anesthetic induced EEG changes (Steyn-Ross et al., 1999; Bojak and Liley, 2005), and (ii) a spatially homogeneous version of the theory has been shown to burst when modified to include a slow modulatory system (Liley and Walsh, 2013).

2. Neural field model for spatiotemporal burst suppression

Here we detail how a neural field model (Liley et al., 2002), subsequently extended to account for the dynamical genesis of the resting EEG and its modulation by anesthesia (Bojak et al., 2004; Bojak and Liley, 2005; Liley and Bojak, 2005; Frascoli et al., 2011; Liley et al., 2011; Bojak et al., 2013), can be plausibly modified to produce bursting-like behavior (Liley and Walsh, 2013), and thus serve as a basis for understanding the emergence of spatially heterogeneous burst suppression seen in cortex. The main advance in this work is that we combine the realistic modeling of isoflurane effects and the extension to a two-dimensional spatial sheet of Bojak and Liley (2005) with an updated version of the slow modulatory system proposed in Liley and Walsh (2013) in order to obtain spatiotemporal activity predictions.

It is perhaps useful to discuss two fundamental limitations of our current approach in advance. First, we are limiting ourselves here to a two-dimensional (toroidal) cortical sheet and use “background” (isotropic and homogeneous) connectivity in order to use computationally efficient activity propagation with partial differential equations (PDEs). Currently there exist a range of mesoscopic approaches available that can incorporate more realistic cortical geometry as well as including “specific” (anisotropic and sparse) connectivity, see for example (Bojak et al., 2010, 2011; Deco et al., 2011; Bojak and Breakspear, 2013; Sanz Leon et al., 2013) and references therein. These typically involve constructing meshes of neural masses and tracking their information exchanges individually. However, such approaches are computationally about an order of magnitude more expensive. Furthermore, to effectively display spatiotemporal pattern formation on a realistically folded cortex is a graphical challenge. We will show here that spatial differentiation of activity emerges even if one uses a simple toroidal cortical geometry with isotropic and homogeneous connectivity. The additional complexity introduced by anatomical folding structures and patchy connectivity are expected to break up long-range coherence further, but should not qualitatively change our results more locally (where the “background” connectivity is a good approximation) and between well-connected but separated regions (where we would expect emergent differentiation). In addition, we use here the well-known “damped wave” propagation PDEs that have been the mainstay of the field since their initial introduction by Jirsa and Haken (1996). It is by now known that one can use “dispersive” propagation PDEs that are more faithful to the actual distribution of axonal fiber velocities (Bojak and Liley, 2010). However, we are using here a parameter set of Bojak and Liley (2005) that delivers realistic EEG activity under the assumption of “damped wave” propagation. Furthermore, the better “dispersive” propagation is also computationally considerably more expensive and technically difficult to implement. Finally, one of the key results of Bojak and Liley (2010) was that more realistic “dispersive” propagation lead to easier spatiotemporal pattern formation. Thus, we expect that the results here would carry over qualitatively to more realistic propagation models, likely showing spatial differentiation earlier on in the burst phase. In summary, we will show here with the computationally simplest model that spatial differentiation in the burst phase can emerge in qualitative agreement with the experimental observation, and we expect that even more realistic modeling will only enhance these emergent effects.

2.1. The (extended) liley model

The electrocortical model of Liley et al. (2002) is constructed at the scale of the cortical macrocolumn. Within each macrocolumn, and extending across all cortical layers, distributed populations of excitatory and inhibitory neurons interact with each other by all possible feedforward and feedback intracortical (local) axo-dendritic connections. Macrocolumns then interact with each other by the exclusively excitatory cortico-cortical (long-range) axonal fibers. The topological organization of this model is well-known, and depicted in Figure 1. In this model cortical activity is described by the spatiotemporal evolution of the mean excitatory h(, t) and inhibitory h(, t) soma membrane potentials. The connection with electrophysiological measurement is through h, which is assumed to be linearly related to the EEG, cf. Bojak and Breakspear (2013). Excitatory and inhibitory neuronal populations are modeled as spatially averaged conductance-based neurons:

Figure 1

Topology of the Liley model (Liley et al., . Its two distinct neural populations (E = excitatory, I = inhibitory) are shown for two separate positions on the cortical sheet. Each one can be considered as representing a single macrocolumn. All synaptic connections that occur in the model are shown by red (excitatory) and blue (inhibitory) disks, respectively. Extracortical inputs to the cortical populations are shown by green fibers. Symbols illustrate the various inputs to the excitatory population in the left macrocolumn and to the inhibitory population in the right macrocolumn, respectively, according to Equations (3, 4).

where ∈ ℝ2 is position on the cortical sheet, subscripts l, k ∈ {e, i} indicate excitatory and inhibitory subpopulations, respectively, and double subscripts represent first the pre-synaptic source and then the post-synaptic target. The parameters hr are the mean resting membrane potentials to which the h decay exponentially with characteristic time scales τ in the absence of inputs I. The fraction in front of the I weighs these inputs, so that the depolarizing effect of additional excitation diminishes linearly and then even becomes hyperpolarizing past the reversal potentials heq, and similarly for the hyperpolarization due to inhibition depending on heq. The weight at the resting potentials is +1 for excitatory and −1 for inhibitory inputs, respectively.

The dynamics of the post-synaptic potentials (PSPs) I are described by critically damped oscillators driven by the mean rate of incoming excitatory or inhibitory axonal pulses A, originally defined as follows (Liley et al., 2002):

For excitatory post-synaptic conductances there are three sources of axonal pulses—local (S), cortico-cortical (Φ), and subcortical (p)—whereas for inhibitory post-synaptic conductances the only source of axonal pulses is local (S), because thalamic and cortical inhibitory axons are essentially short-range on the basis of existing neuroanatomical evidence. For these equations, at a given location a single pre-synaptic (Dirac delta) spike A(t) = δ(t) would produce a so-called “alpha function” response

with the Heaviside step function Θ(t). The alpha function α is normed to one for integration over time, hence the pre-factor in Equation (5) is proportional to the charge transfer of the induced PSP. Furthermore, δ = 1/γ is the characteristic time scale of the PSP's exponential decay. Since I(t = δ) = Γ is the maximum amplitude of the PSP, δ is also the rise time to peak amplitude. Note that since we have collapsed all cortical layers into one sheet without radial extension, this rise time is that of the PSP conducted to the soma rather than at the synapse in the dendritic tree. Conduction through a passive dendritic cable effectively leads to a “flattened” PSP at the soma with lower maximum amplitude and prolonged rise and decay times.

The S and S are respectively the local mean excitatory and inhibitory firing rates, and are assumed to be instantaneous sigmoidal functions of the h of the form

and the Nα and Nβ factors in the A above multiply these local firing rates by the number of synaptic connections formed with the target populations. The Smax are the maximum mean firing rates, and the μ and σ can be understood as the mean and standard deviation, respectively, of the firing thresholds of the populations, which are taken to be roughly normally distributed. The propagation of axonal pulses by the excitatory cortico-cortical fiber system Φ is described here by the following well-known “damped wave” equation (Jirsa and Haken, 1996; Robinson et al., 1997; Liley et al., 2002; Bojak and Liley, 2010):

But for the λ terms, this would be an inhomogeneous wave equation with conduction velocity v, propagating the local excitatory firing rate S. However, due to these terms the wave gets suppressed roughly exponentially with distance with a characteristic spatial scale λ.

Finally, there is also extracortical synaptic input in the form of the p. These inputs can be considered to be mainly due to thalamic afferents. If the p were constant, then for the model parameters chosen here the system would quickly converge to a static equilibrium point. It is hence the imposition of noise on these inputs which effectively drives the neural activity. This noise is taken to represent the average over the varied extracortical synaptic input to the many thousands of neurons in a neural population, for the case in which there is no strong external (sensory) drive that would lead to clear correlations of the synaptic inputs to the individual neurons. We follow here essentially the approach of Bojak and Liley (2005) for noise generation. Thus, for the sake of computational simplicity noise is imposed only on p, whereas p is taken to be constant. At every grid point of the two-dimensional cortical sheet normally distributed noise is generated independently, but with the same mean p, and a standard deviation that is 10% of this mean. However, we filter this noise spatiotemporally, both to achieve more biological realism and to make it easier to achieve numerical stability. We follow the Fourier space procedure of Bojak and Liley (2005) for the spatial filtering, but use the Catmull-Rom spline procedure detailed in Bojak et al. (2011) for the temporal filtering, with lowpass −3 dB points at 75 Hz and 2/cm, respectively (Bojak and Liley, 2005). Thus, the noisy input oscillates equally at all frequencies, but only up to about 75 Hz, and is identical for neighboring grid points, but becomes uncorrelated at cortical distances greater than about 0.5 cm. The spatiotemporal p noise breaks the otherwise perfect homogeneity and isotropy of the system, and consequently acts as seed for the heterogeneities observed in the burst suppression phase. However, the characteristics of the spatiotemporal structures that emerge in the burst suppression phase do not otherwise depend on the noise; hence in particular they do not depend on the details of the noise filtering, and can be elicited with white noise driving.

The model we have described so far has to be extended for a realistic description of the effect of general anesthetic action. In particular, the effect of isoflurane on the rise time δ of the PSPs from zero to maximum amplitude Γ, and on the subsequent decay time ζ back to Γ/e (measured here from the start of the PSP, not from the peak) can be parameterized as follows in the form of a Hill equation (Bojak and Liley, 2005):

where c is the aqueous concentration in mM. Thus, the main effect is a prolongation of the decay of the inhibitory PSPs. In addition, the maximum amplitudes of the PSPs also diminish with increased isoflurane concentration, which is also the case for excitatory PSPs:

While for consistency with laboratory based estimates we choose to parameterize isoflurane level in terms of its aqueous concentration, it is important to appreciate that because isoflurane is a volatile gas, clinically its level is typically reported in terms of its concentration in the expired air (which is assumed to be in equilibrium with the blood and hence the extracellular fluid of cortical neurons). At normal body temperature 1.3% isoflurane are equivalent to an aqueous concentration of about c ≃ 0.27 mM (Franks and Lieb, 1996). Typical isoflurane concentrations encountered in clinical anesthetic practice range from 0 to 2% of the expired air, equivalent to aqueous concentrations 0 − 0.42 mM (Mapleson, 1996). A measure commonly employed in anesthetic practice is the minimum alveolar concentration or MAC of an anesthetic agent. It is essentially defined as the concentration of gas in the lungs required to prevent movement in 50% of subjects in response to a painful surgical stimulus. In the case of isoflurane 1 MAC ≃ 1.17% ≃ 0.243 mM for an adult at normal body temperature (Mapleson, 1996). It should be noted that in combination with other anesthetic agents like nitrous oxide, less than 1 MAC of isoflurane will be required to reach the 50% end-point.

It is straightforward to introduce Equations (12, 13) to the Liley model by changing the Γ according to isoflurane concentration, i.e., Γ → Γ(c) = Γ0H(c) with the Γ0 now having the same values as the Γ had in the standard Liley model. However, Equations (9–11) are more problematic. The decay time of the alpha function in Equation (5) changes linearly with its rise time, thus one cannot match the experimental result that only the decay time is prolonged under anesthesia. Consequently, the following modification of Equation (2) was introduced (Bojak and Liley, 2005)

Notably, for ε → 0 one finds that , and γ → 1/δ with a removable discontinuity. Defining these variables as continuous with the limit, the new Equation (13) then becomes identical with the old Equation (2) in this limit. The corresponding response to a single pre-synaptic spike A(t) = δ(t) now becomes a bi-exponential function

where we have suppressed the concentration dependence. Again the pre-factor in Equation (16) is proportional to the charge transferred, since β(t) is normed to one for integration over time. Note that now I (t = δ) = Γ, which for ε = 0 becomes the previous result since in this limit again δ = 1/γ. More generally, for ε → 0 we have β(t) → α(t) at all times and the alpha function is the “sharpest” response β(t) ≥ α(t). Clearly with this new form we can keep the rise time parameter δ constant, while changing the ε so as to achieve a desired decay time ζ. Given the changes imposed by isoflurane in Equations (9–11), one can solve for the appropriate ε numerically in dependence on the concentration c. However, here we will use the excellent approximation formula presented in Liley et al. (2011), which can be written as

where W−1 is the −1 branch of the Lambert-W function. Here we assume that κ = 1, thus ε = 0, and only the inhibitory decay time is affected.

Equations (1–4, 7–8) represent a system of eight coupled non-linear PDEs that define the standard Liley model. Changing the PSPs of Equation (2) to those of Equations (14, 15) defines the extended Liley model. It is therein understood that 1/δ of the extended Liley model equals γ of the standard one, so that for ε → 0 both become identical. Finally, Equations (9–13) parameterize the effect of isoflurane on the extended Liley system. Here Equations (12, 13) can be used straightforwardly as determining Γ(c), but in order to use Equations (10–11) one additionally needs Equation (18) to translate them into changes of the γ(c) and parameters. In Bojak and Liley (2005) extensive parameter searches were performed. All selected parameter sets gave rise to a plausible resting EEG power spectrum (‘1/f’ low frequency activity with an alpha peak in the 8–13 Hz range) under noise driving, retained a stable equilibrium point for increasing isoflurane concentration and hence remained in a quasi-linear dynamical regime, and showed the experimentally observed drop of the alpha peak to low frequencies for increasing isoflurane concentration. Some parameter sets furthermore exhibited a so-called “bi-phasic” transient surge in total power during simulated anesthesia induction, as observed in several experiments (Kuizenga et al., 1998, 2001). The parameter values used in this paper correspond to one of these “bi-phasic” parameter sets, and are listed in Table 1.

Table 1

Mean population parameter values used to obtain bursting in the Liley model.

Definition	Excitatory (target)		Inhibitory (target)
Passive membrane decay times	τe	65.815 ms	τi	130.13 ms
Resting membrane potentials	hre%	−78.422 mV	hri%	−72.959 mV
Maximum firing rates	Smaxe	0.39535/ms	Smaxi	0.15439/ms
Firing thresholds (FTs)	μe	−51.656 mV	μi	−47.267 mV
Standard deviations of FTs	σe	2.8669 mV	σi	4.3250 mV
Synaptic recovery times	τrece	800.00 ms	τreci	600.00 ms
Synaptic depletion factor	fe	1.2500	fi	0.17500
EXCITATORY SOURCE
Reversal potentials	heqee	−5.7891 mV	heqei	−1.6566 mV
PSP peak amplitudes	Γee	0.18424 mV	Γei	1.8771 mV
PSP rise times to peak	δee	9.1059 ms	δei	1.2103 ms
Number of intracortical synapses	Nβee	3410.8	Nβei	2738.9
Number of cortico-cortical synapses	Nαee	3616.3	Nαei	2905.1
Cortico-cortical decay scale	λee	24.000 mm	λei	24.000 mm
Cortico-cortical conduction velocity	vee	2.1042 mm/ms	vei	2.1042 mm/ms
Rate of extracortical input	pee	9.3193/ms	pei	3.1563/ms
INHIBITORY SOURCE
Reversal potentials	heqie	−86.675 mV	heqii	−84.596 mV
PSP peak amplitudes	Γie	1.5969 mV	Γii	1.0838 mV
PSP rise times to peak	δie	2.5985 ms	δii	9.6946 ms
Number of intracortical synapses	Nβie	863.89	Nβii	267.92

2.2. Slow and activity-dependent synaptic bursting mechanism

In this work we consider receptor desensitization and synaptic vesicle depletion during periods of high neuronal population activity, and the homeostatic recovery of synaptic readiness during periods of low neuronal activity, as the slow mechanism that can modulate the excitability of cortical tissue. Such activity-dependence of synaptic efficacy has been observed during burst suppression induced with isoflurane (Kroeger and Amzica, 2007). In practice, we will modify the maximum PSP amplitudes Γ that can be obtained, which directly depend on the available pre-synaptic amount and post-synaptic impact of neurotransmitter. Instead of considering these quantities as parameters as in the extended Liley model (where they act as control parameters that can be changed according to the concentration of an anesthetic agent), we now consider them as variables with their own slow dynamics coupled to the neural activity. Our ansatz is a common phenomenological model for activity-dependent synaptic depression (Bressloff, 2012). It represents a rate-based version of the model proposed by Tsodyks and Markram (1997), under the assumption that the processes responsible for the recovery of synaptic efficacy evolve on a time scale much slower than those associated with that of synaptic depletion (e.g., receptor desensitization and synaptic vesicle depletion):

with l, k ∈ {e, i} indicating the excitatory and inhibitory subpopulations and S(h) is the local population firing rate of Equation (7), as before. Local neurotransmitter depletion is here considered to be directly proportional both to the strength of the PSPs, represented by Γ itself, and to their frequency, represented by S. In the absence of neural activity S = 0/s, there will be an exponential return of Γ to the resting value Γr with a characteristic recovery time τrec. However, if there is no recovery τrec → ∞ and we have constant neural activity S > 0/s, then Γ will exponentially decay to zero with a characteristic depletion time 1/(ρdep S). Note that we have assumed that in the pre-synaptic recovery and decay there is no dependence on the target (on the index k), since these processes will be determined by the activity of the source. However, in the post-synaptic impact on the maximum amplitude of the PSP, we allow a dependence on the target, since the response will depend on the morphology and physiology of the receiving neurons.

Now consider the case where we have both depletion and recovery. We will choose some homogeneous h(, t) ≡ h0 so that S(h0) > 0/s. The synaptic system will then converge everywhere to an equilibrium value easily calculated by setting the left hand side of Equation (19) to zero:

This equilibrium value is always smaller than the resting one, i.e., f > 0. We can scale our ansatz in terms of this equilibrium value and then obtain

This scaling conveniently removes the post-synaptic dependence from the dynamical equations. Hence if we assume that the scaled initial conditions are identical Γ(, t = 0)/Γ0 = Γ(, t = 0)/Γ0, then in practice we only have to solve the two equations of Equation (23) instead of the four of Equation (19), obtaining the four potentially different peak amplitudes via the scaling in Equation (22). Our choice for the initial conditions is to start them all at equilibrium Γ(, t = 0) = Γ0, thus C(, t = 0) = C(, t = 0) = 1, and then we can use these reduced computations with subsequent scaling. The only constraint we have imposed on the mean membrane potentials h0 here is that they should lead to non-zero population firing rates, which however is always the case unless one assumes unphysiological infinite polarization. The depletion coupling constants between the neural activity and the peak amplitude ρdep are empirically unknown, and would have to be determined laboriously from observations of dynamical changes of synaptic efficacy. However, in terms of the model proposed here, if one specifies the equilibrium values of the soma membrane potentials h0 and how much the neurotransmitter reservoir is depleted at the corresponding activity levels (Γ0 vs. Γr), then this in turn determines the depletion coupling constants

In practice we make an implicit choice of the coupling constants by choosing the f for the system.

Now we wish to combine this synaptic system with the extended Liley model for anesthesia. For the parameter sets provided by Bojak and Liley (2005), that model has stable equilibrium points. That is to say, if we re-write the extended Liley model in the abstract form

with a suitable differential operator , a function F and a noise drive P, then there exists a solution

so that for P(, t) = 0 the system is static. Furthermore, since this equilibrium is stable, after small and transient disturbances the system will return dynamically to *. We now make the following replacements in the extended Liley system

i.e., we replace the parameter values Γ of the extended Liley model with the variables Γ(, t) of the synaptic system that we have just described. Together with the coupling to the neural activity explicit in Equation (23) this closes the combined system, which we will call the bursting Liley model henceforth. We now make the following convenient choices

This homogeneous initial state of the bursting Liley model must now be an equilibrium point by construction: while Equation (27) is calculated with the parameters Γ, we have arranged it so that the equilibrium value Γ0 of the synaptic system at the resulting mean soma membrane potentials has the same value as that parameter. This is simply achieved by fixing the Γr for a given f according to Equation (20), i.e., Γr = Γ0(1 + f). Hence the equilibrium of one system is compatible with that of the other, and if we start them off in their respective equilibrium states nothing will change. However, there is no guarantee that this constructed equilibrium point of the bursting Liley model will be stable.

Previously, we had incorporated the effects of isoflurane into the extended Liley model in part by replacing the standard parameter Γ0 according to Equations (12, 13) with the anesthesia-dependent Γ(c) = Γ0 H(c). In the bursting Liley model these parameters have become state variables with their own dynamics due to synaptic depletion and recovery. Hence the synaptic dynamics pertaining to the (pre-synaptic) source likewise must be multiplied by the anesthesia-dependent H(c) to compute the (post-synaptic) amplitude induced at the target:

However, since the synaptic dynamics are now coupled to the spatially variable cortical activity, we need to adjust our synaptic inputs to Equation (14):

Here the first term in A is multiplied with C(, t) at the same time and position, since it represents local and quasi-instantaneous synaptic input. For the second term of A, the C(, t) term is instead included through Equation (33). The right hand side of this propagation equation, while written in terms of the (, t), effectively encodes the signal at a distance location ′, sampled there at time t′, and then transported with velocity v to the local position with a conduction delay t − t′; see for example (Bojak and Liley, 2010) for an explanation in terms of Green's functions. Thus, we now propagate the firing rate as scaled by the pre-synaptic efficacy of the neural populations at a distant position ′ at the time t′. We note that this is not quite physiologically accurate either, since the synaptic dynamics should be evaluated at (, t), not (′, t′), albeit driven with the firing rates from ′ delayed by t − t′. This could be achieved by setting NβC → NβC and NαΦ → Nα CΦΦ in Equation (31), removing C in Equation (33), and then tracking separately the spatiotemporal dynamics of C and CΦ, respectively, where the latter would have Φ instead of S in Equation (23). However, this would double the effort for computing the synaptic dynamics and could have potentially undesirable consequences for the separation of distant sources, see the Discussion for further detail. Finally, for simplicity we have assumed here that the extracortical input p remains unchanged. This assumption will likely need to be improved upon for greater physiological realism, i.e., we do expect that in particular thalamic activity will also be modified by anesthesia. However, our current focus on only the cortical side allows us to highlight the proposed bursting mechanism without potential interference from complex interactions between extracortical and cortical structures. Given the assumption that p is constant (or in the case of p, that its mean is constant), we would not expect the pre-synaptic efficacy to change. Overall, if we switch off the synaptic dynamics C ≡ 1, we recover exactly the extended Liley model of Bojak and Liley (2005).

The bursting Liley model hence consists of Equations (1, 7, 14, 15, 21–23, 30–33), with the influence of isoflurane anesthesia being parameterized by Equations (9–13, 18). In practice we choose the f > 0, and as before use the “combined fixed point” initial state of Equation (29). The parameter values we have chosen are listed in Table 1. For the synaptic system, we have followed qualitatively the work of Tsodyks and Markram (1997) in assuming a possible range of about 250 ms to 1000 ms for τrec, and values between 0.1 and 2.0 for f. The values used in this paper were chosen after computational experimentation with various settings, and were selected because they lead to bursting only for relatively large concentrations of isoflurane. Clearly, more systematic and comprehensive scans of the available parameter space and better understanding of the dependence of the observed dynamics on these parameter values are needed in order to elucidate the mechanisms proposed here. However, it takes considerable computing time to simulate such large spatial systems. In order to accomplish a proper analysis of the parametric dependencies, one will likely need to find approximate but rapid evaluation methods, similar to replacing the full simulation with an eigenvalue calculation as in Bojak and Liley (2005). The development of such methods is beyond the scope of this article, here we want to demonstrate in a pilot study that we can qualitatively reproduce the spatial differentiation in burst suppression that has been observed experimentally.

2.3. Numerical simulations

All our simulations are performed on a two-dimensional cortical sheet discretized by a 512 × 512 numerical grid, where we assume a grid spacing of Δx = Δy = 1 mm. The effective simulation area of 2,621.44 cm2 corresponds roughly to the size of an entire human cortex (Im et al., 2008). Smaller grids, with or without larger grid spacing, have been used to investigate parameter dependencies more rapidly, but the results presented in this paper were all obtained on this standard grid. In order to avoid boundary effects we have made the numerical grid toroidal, i.e., if we number the grid points 0–511 along one dimension from left to right, then the grid point to the left of 0 is 511, and the grid point to the right of 511 is 0, and this is true for both dimensions. Obviously such a geometry is artificial as compared to the real brain. However, since it leaves all numerical grid points entirely equivalent, this together with the isotropic and homogeneous “background” connectivity implicit in the PDE propagation makes minimal assumptions about the actual geometry and specific connectivity of the brain. Basically, it represents a kind of anatomical “null hypothesis” from which any anatomical detail will deviate; and the more homogeneous the brain turns out to be in an effective sense, the better this approximation will represent its activity. As argued above, at significantly increased computational costs one can improve this description with neural mass meshes, but this is not expected to change the results obtained here at least qualitatively.

The only dependence on space is found in the propagation PDE of Equation (33). Hence the other dynamics are effectively described by a set of independent ordinary differential equations (ODEs) in time at every grid point. We solve all these ODEs with the following simple method: First, any higher time derivatives are turned into first derivatives by defining auxiliary variables, e.g., d2g/dt2 = f(g) becomes and . Next, we solve these first order ODE systems with the forward Euler method. Obviously many more efficient numerical schemes exist. But in our experience they occasionally fail for specific parameter settings with the Liley model, whereas the forward Euler method always remains stable. Thus, we trade speed for guaranteed stability here. In the propagation PDE, the Laplacian is approximated by a five point stencil, i.e., to estimate the Laplacian at a grid point, we use the value at that point and those of its four horizontal and vertical nearest grid neighbors. We find that numerical stability is increased, if in this PDE we likewise estimate the second derivative in time directly by considering the current, previous and future values (and solving for the future one), rather than first rewriting them into first order derivatives as for the other dynamics. We use MPI-C to parallelize the computation across multiple nodes (threads and/or cores). This involves splitting up the grid into patches assigned to the individual nodes. We note that since the only spatial dependence in the dynamics arises from the Laplacian, and since we approximate it with a five point stencil, the only required communication between these nodes is that of the proximate part of the one grid point deep boundary of the local patch to the nodes working on the adjacent patches. This limited need for communication between nodes allows for very efficient parallel computation.

How the noise driving the system is generated and filtered in a mathematical sense has been described above, here we will add the following technical comments: The temporal Catmull-Rom spline filter is obviously local in space, and hence in a parallel setting can be performed by every individual compute node on the grid points assigned to it. However, the initial noise generation is done in Fourier space, to allow spatial filtering by a simple multiplication at every (Fourier space) grid point, followed by an inverse Fourier transformation. We use FFTW (Frigo and Johnson, 2005) to perform the inverse Fourier transform in parallel across the available compute nodes. This means that the random number generation and the Fourier space filtering can be done local in each node on its part of the Fourier grid, while FFTW organizes the communications between the nodes involved in the inverse Fourier transform.

The time step used in our simulations is Δt = 5 · 10−5 s, which for our chosen grid spacing is sufficient to achieve stable and convergent results. However, we save the simulation results neither at every time step, nor the entire system state, nor at the internal double floating point precision (8 Bytes). The reason is that the 10 state variables of the bursting Liley model saved for 512 × 512 grid points at 8 bytes per point already would require 20 MB of hard disk space per time step, and thus at full time resolution a mere 3 s of run time would generate more than 1 TB of data. In practice, we typically save between one and four selected state variables with 250 Hz, converting to single floating point precision (4 Bytes) in the output. This still generates data files of many GB for our longest runs. Finally, even this data reduction is not sufficient to generate suitably sized animations of our results. Basically, the fine detail of a 512 × 512 grid leads to low compression efficiency of the employed (H. 264) movie codec. For producing animations we hence tile the output grid into squares of 4 × 4 grid points and average over these to obtain an effective 128 × 128 grid with smoother values that compress better. This explains the mild visual disparity between our figures (at full 512 × 512 resolution) and the animations, even though they are produced from the same underlying data set. We also produce video frames at an even lower sampling rate in time, and we use variable sampling rates to selectively speed up uneventful parts of the video. A time counter in the videos keeps track of the sampling rate, and occasional choppiness and blurriness in the videos does not reflect actual discontinuities in the simulations but merely low sampling rates and aggressive video compression.

3. Results

We explore the influence of isoflurane on the model in a long simulation run presented in Figure 2. The entire simulation also has been animated as Movie 1, included in the Supplementary Material. In Figure 2A we show the time course of the isoflurane concentration that we have imposed. First the system is run free of anesthesia (0 MAC) for 10 s. We call this the first plateau in the following. The equilibrium values of the system are used as initial conditions. Hence there are no transient dynamics, which allows us to estimate a power spectral density (PSD) from the h time series. Then we increase the concentration linearly to 0.5 MAC (equivalent to 0.1215 mM or 0.585% inspired at normal body temperature) over 10 s, and keep the system at this concentration for another 10 s. This second plateau corresponds to a light anesthesia state, without burst suppression, and again we can estimate a PSD here. Next we increase the concentration linearly to 1.0 MAC (equivalent to 0.243 mM or 1.17% inspired), and keep the system there for 40 s. This third plateau corresponds to a state of deep anesthesia, with burst suppression, and we can estimate a PSD here as well. After that, we increase the isoflurane concentration again for 10 s to 1.5 MAC (equivalent to 0.3645 mM or 1.755% inspired), and maintain it at this value for 10 s. Bursting is abolished at this fourth plateau, and we compute another PSD here. Finally, we raise the concentration for another 20 s up to 2.5 MAC (equivalent to 0.6075 mM or 2.925% inspired). This demonstrates that the system has finally returned to a regime without bursting.

Figure 2

(A) Imposed concentration of isoflurane (red curve), and the h response (blue curve) at the cortical location indicated by black arrows in the snapshot panels below. Different plateaus of concentration are labeled “C,” “D,” “E,” and “F.” Arrows point to the central times of the corresponding time series shown below. (B) PSDs of h averaged over the entire grid and normed to unit area for plateaus “C” (blue), “D” (green), “E” (red), and “F” (cyan). The motion of the alpha peak to lower frequencies persists qualitatively into the burst suppression phase “E” at much increased power. (C1) Snapshot of the h activity of the cortical surface at 0 MAC isoflurane. The size of h is indicated by both height and color, cf. the color bar. A black arrow shows the position from which the corresponding time series were recorded. (C2) Time series of h (blue) and Γ (green) over the 10 s of the “C” plateau. Regular alpha rhythms in h and slow Γ oscillations around the standard value Γ0 can be seen. (D1) Snapshot at 0.5 MAC. (D2) Time series of the “D” plateau. The oscillations of h have larger amplitude at a lower average. The slow Γ oscillations now occur at an elevated level. (E1) Snapshot at 1 MAC. Burst suppression patterns have emerged and move across the cortical surface. (E2) Time series of the “E” plateau. Burst suppression is apparent both in h and Γ, with a rapid drop in Γ caused by the strongest h oscillations. An animation of this simulation is provided as Movie 1 in the Supplementary Material.

In Figure 2B, we see PSDs estimated over these plateaus. At each concentration, we have calculated PSDs for every individual grid point from the h time series for the entire duration of the plateau, using a Welch estimate with a 2.5 s window and 50% overlap, and then have averaged the resulting 262,144 PSDs. We have normed these average PSDs to have unit area, i.e., their total power over all frequencies is one. This makes it easier to compare them visually. Please note that the parameter set used here is a “bi-phasic” one, see the discussion in Bojak and Liley (2005). Consequently, the total power at the second plateau is actually increased over that at the first plateau by a factor of 1.26. We see the characteristic shift of the alpha resonance to lower frequencies, in this case initially accompanied by a sharpening of the peak. Without the slow synaptic system as in Bojak and Liley (2005), a further increase of the isoflurane concentration would move the former alpha peak to ever lower frequencies, accompanied eventually by strong damping of the peak and a reduction of the total power, completing the bi-phasic power change. However, with the introduction of the slow synaptic system we see a burst suppression pattern emerge, and the large amplitude oscillations imply a drastic increase in total power by a factor 131 over the resting values. Yet we see that the PSD obtained from the third plateau—as far as frequency content is concerned—roughly follows what is expected without the synaptic system: The majority of the power, which is generated by the large bursts, is located where one would expect to see the former alpha resonance in the previous model of Bojak and Liley (2005). In other words, the bursts roughly conserve the regular dynamics of the system, in particular of the former alpha resonance. At the fourth plateau the system has ceased to burst but still shows elevated total power 1.44 times larger than at rest.

In Figure 2C1 we see the system at rest under noise drive. The visible structure is hence basically that of the spatial correlations we have included in the noise. The corresponding time series in Figure 2C2 shows the typical waxing and waning of a resting alpha rhythm in h. Unsurprisingly, Γ oscillates slowly at values about the equilibrium value of the slow synaptic system Γ0, cf. Table 1. In Figures 2D1,D2 we see the corresponding state at light anesthesia. The overall h is now lower across the grid, indicating smaller firing rates on average. However, as we can see in the time series the amplitude of the oscillations has increased, corresponding to the power increase expected for this “bi-phasic” parameter set. The oscillation frequency also has become lower, though this is easier to see in the PSDs of Figure 2B. We see that Γ is still oscillating slowly, but around values somewhat higher than Γ0, because the synaptic resources are not as rapidly depleted by the reduced excitatory firing rate of the depressed h.

In Figure 2E1 burst suppression patterns have emerged. These patterns are clearly independent of the noise drive. In the Movie 1, included in the Supplementary Material, one can see how this is a snapshot of “burst waves” moving across cortex, with centers of burst activity spontaneously forming and disappearing. The geometry of these excitations is complex and constantly changing. We see in Figure 2E2 that the strongest oscillations in h are associated with a rapid drop in Γ due to the synaptic depletion during these periods of high firing. This lowering of Γ quickly suppresses the burst by reducing the self-excitation of cortex. This is then followed by a recovery to values of Γ that are large compared to those at rest or light anesthesia. This recovery to high values of Γ is driven by the isoflurane-induced reduction in the mean excitatory firing rate during the suppressed periods. In turn, these strong PSPs eventually destabilize the neural system, leading to another burst. The burst suppression pattern persists as the anesthesia concentration is being increased again, up to quite high concentrations. In Movie 1, one can see that the amplitudes of the “burst waves” eventually become smaller and smaller until the burst activity fades away into regular noise driven activity. We do not show here corresponding plots for the fourth plateau with abolished bursting, but they would look similar to Figures 2D1,D2 with a further reduced mean value h, even lower oscillation frequency, and a Γ that is on average even higher.

The burst suppression phase shown in Movie 1 of the Supplementary Material makes obvious that one cannot expect global synchrony of the burst suppression across cortex. A multitude of transient spatiotemporal patterns emerge, travel across cortex, and dissolve. This is also shown in Figure 3, which shows time series from three well-separated locations on the simulated cortex at a specific point in time. There is little evidence of strong systematic correlations. While one might expect that the propagation of “burst waves” should lead to correlations with temporal delay at these distances, other burst features emerge across these spatial scales and interfere with the burst timing. Without observing spatiotemporal pattern globally, it hence will be difficult to find systematic correlations of the bursts at large distances. However, locally it may be possible to track the regular motion of burst patterns, e.g., at a point close to the one labeled “a” one might see bursting appear with a delay, characteristic for the “burst wave” passing through these two points sequentially. Overall, we expect stronger synchronization—or at least consistent phase differences from traveling patterns—at shorter distances, whereas at longer distances such correlations will be basically accidental. Thus, one would expect to see considerable spatial differentiation if one records from several spatial locations, as in Lewis et al. (2013). How many electrodes would be seen to burst at the same time would depend on the size and motion of the emerging spatiotemporal burst patterns.

Figure 3

(A) Snapshot of the h activity of the cortical surface at simulation time 54.94 s under the influence of 1 MAC isoflurane. Black arrows with labels “a,” “b,” and “c” point to the cortical locations of the time series shown in the other panel. Note the toroidal boundaries, e.g., the circular burst front that appears cut off around (x,y) = (44.8, 0) cm continues at x = (44.8, 51.2) cm. (B) Time series of h (blue) and Γ (green) taken from the three different positions marked as “a,” “b,” and “c” in the other panel. It is obvious that bursts are not generally synchronized in time at the different positions. Spatiotemporal correlations from propagating “burst waves” can occur, but are removed at larger distances by the interference from other emergent patterns. An animation of this simulation is provided as Movie 1 in the Supplementary Material.

Local variation of cortical tissue properties, reflected in the model evaluation by a change in the parameters, may also affect the ability of some part of cortex to participate (fully) in the spatiotemporal burst suppression dynamics. Such variation of tissue properties can be natural and develop intrinsically, or could be induced extrinsically by physical insult or the application of drugs. We have seen that bursts are associated with slow but large oscillations in the excitatory peak amplitudes of the PSPs. It is important to note that there are two different effects determining the general size of the Γ. On one hand, anesthesia is reducing Γ directly as parameterized by the Hill factor H(c), cf. Equation (30). On the other hand, the reduction in the average h, mostly due to the strong prolongation of the inhibitory PSPs with anesthesia, means that the average excitatory firing rate S decreases. This in turn leads to less synaptic depletion and hence actually a rise in Γ, cf. Equation (19). The net effect with increasing concentration is actually an increase of Γ, and this is crucial for the onset of bursting. If one increases anesthesia further, eventually the Hill factor begins to dominate and Γ decreases again.

The same can be said for the Γ, and the corresponding balance between the Hill factor H(c) and the reduction in S for increasing anesthesia. However, we see that in the standard parameters the excitatory synaptic depletion factor f = 1.25 is much larger than the inhibitory one f = 0.175. This means that there is much less room for Γ to grow, since the steady maximum is Γr(c) = Γ0(1 + f)H(c). One simple idea for reducing the ability of cortical tissue to participate in burst suppression is hence to increase the growth of inhibition with anesthesia by raising f. What do we expect to be the effect of this increased inhibition, in particular concerning the excitatory Γ? In general we expect h and h to decrease even more rapidly with increasing concentration of anesthesia, due to the boosted inhibition. But silencing the cortex also decreases synaptic depletion, so we actually expect a stronger initial Γ growth for increasing anesthesia. It is hence not a priori clear whether the more rapid decrease in h or the more rapid increase in Γ dominates, and in consequence whether bursting is abolished or maintained, respectively. We note that for our regular parameter set bursting is abolished at higher concentrations even though Γ is still increasing, because h has then decreased too much. We may hence expect that an increase of f alone can stop the bursting.

To test this, we use a typical simulation at 0.25 mM isoflurane. However, in a circular patch of tissue we set f = 1.25 instead of the standard f = 0.175, while leaving this parameter at f = 0.175 across the rest of the cortical sheet. Using Γr(c) = Γ0(1 + f)H(c), we have inside this patch Γr(0.25) = (3.5174, 2.3872) mV and outside Γr(0.25) = (1.8369, 1.2467) mV, respectively, while everywhere Γr(0.25) = (0.37703, 3.8414) mV. As shown in Figure 4, outside of the circular patch burst suppression patterns emerge as usual, see Figures 4A,B, while in the dead center of the circular patch there is no sign of such activity, see Figures 4A,D. Hence inside the patch the greater decrease in h was more effective than the greater increase in Γ, compare Figures 4B,D. Interestingly, at the border of the circular patch, see Figures 4A,C, we see largely the same state as for the center, but there appear to be some “quasi-bursts”. Actually, this is activity spilling into the circular patch from the outside through the propagation with Equation (33). The characteristic spatial decay scale of this propagation is λ = 2.4 cm. Given a radius of 9.6 cm of the circular patch, we expect a signal from the outside to have fallen to less than 2% of its original value at the center. So it is unsurprising that any outside influence on the center is not obvious to the eye, but that close to the rim we see stronger echoes of the surrounding burst activity. In Movie 2 in the Supplementary Material h (top panel) and Γ (bottom panel) animations are shown. Here one can observe the bursting waves collide with the circular patch, and then fade as they penetrate deeper. We note that an increase to for example f = 0.5 in the patch is not sufficient to abolish bursting in this manner, illustrating that it is the balance between the decrease in h and the increase in Γ which determines whether self-excitation is possible.

Figure 4

(A) Snapshot of the h activity of the cortical surface at simulation time 16.78 s under the influence of 0.25 mM isoflurane, where in a circular patch (center (x, y) = (20.0, 20.0) cm, radius 9.6 cm) the inhibitory synaptic depletion factor f has been increased from 0.175 to 1.25, leading to an increase of Γr there by a factor 1.91. The size of h is indicated by both height and color, cf. the color bar. One can see that the circular patch does not participate in the burst suppression pattern. Black arrows with labels “B,” “C,” and “D” point to the cortical locations used in the other panels. (B) Six seconds long time series of h (blue), h (green), Γ/Γ0 (red), and Γ/Γ0 (cyan) around the time of the snapshot from a point outside of the circular patch. Burst suppression is clearly visible in all variables. (C) Time series from just inside the circular patch. There is no local burst suppression, but some of the outside burst activity spills in. (D) Time series taken from the center of the circular patch. There is neither local burst suppression nor spill-in. Animations of this simulation are provided as Movie 2 in the Supplementary Material.

The dynamics of Γ are of course much slower than those of h, and we can track the h burst fronts by the progression of the lowest dips and valleys in Γ. This corresponds to rapid synaptic depletion in high firing regions. The extent of “spill-in” from the outside into the circular patch is also easier to discern in Γ: the center of the patch remains at roughly constant values, while at the rim Γ drops when high firing propagates into the patch. We note that the strong oscillations that one can see as Γ drops rapidly in Figure 4C show up in the movie as a kind of “bouncing” rather than a smooth advance of the burst fronts. To understand this better, we provide Figure 5. It shows part of the time series labeled “a” in Figure 2B. However, in addition to h (blue curve) and Γ (green curve), it also shows the mean excitatory firing rate S(h) as red curve. So that one ordinate can be used for both Γ and S(h), we have normalized the latter by the maximum excitatory firing rate S(h)/Smax. The basically symmetric oscillations of the mean membrane potential h around an average value translate into strong “spikes” in the mean firing rate S(h). This is due to the sigmoidal nature of Equation (7), combined with the fact that the average h is about 5.5 standard deviations σ = 2.8669 mV below the average firing threshold μ = −51.656 mV, leading to low firing rates. Only the strongest depolarizations in the burst come close to this threshold—though even they do not quite reach it here, as we can see, since S(μ)/Smax = 0.5 by definition. Thus, the relationship between h (local field potentials and EEG) with firing rates is highly non-linear in the anesthetic regime. It is obvious from Figure 5 that the jagged drop of Γ is simply caused by strong synaptic depletion induced by “spikes” in the mean excitatory firing rate.

Figure 5

Time series of the mean excitatory soma membrane potential . Note that both Γ and S(h)/Smax map by value to the black ordinate on the right, though with different units. The time series shown here is part of the times series labeled “a” in Figure 3B. One sees that the strongly non-linear relationship between h and S in the anesthetic regime transforms the “symmetric” h oscillation that would be visible in local field potentials and the EEG during the burst phase into strong “spikes” in the firing rate S, and consequently to a “jagged” appearance of the synaptic depletion of Γ.

Finally, we also considered the influence of the spatial scale of brain connectivity on the spatiotemporal expression of burst suppression. As mentioned above, in this simplified model it is represented by the parameter λ, the characteristic length scale of the exponential decay of activity propagated with Equation (33). Its regular value according to Table 1 is λ = λ2 = 2.4 cm. This is a length scale one might associate with a brain region and cortico-cortical connections, in particular since the influence of activity at a point would be felt across a distance of several λ. We vary this length scale up λ = λ1 = 2.7 cm and down λ = λ3 = 2.1 cm to investigate the impact of brain connectivity on the dynamics. In Figure 6 we see the dependence of the spatiotemporal activity on adjusting this parameter. The spatial extent of the emerging burst patterns clearly becomes smaller as λ is being decreased, cf. Figures 6A–C. This is particularly obvious in the corresponding animations for the different λ values, Movies 3–5 of the Supplementary Material, where we can see that at λ1 large parts of cortex are recruited in the bursts, whereas for λ3 bursting is much more localized. Our standard λ2 represents an intermediate case. However, the time interval between bursts for these different λ appears at first sight comparable, see Figure 6D. To be more quantitative, one can use the point where Γ drops lowest as a convenient marker for the time of a burst peak. To carry out automatic computations for 50 s time series for every grid point, we select the deepest minimum within a specific continuous “burst peak region” defined by Γ ≤ 0.05 as burst peak time, and we remove inter-burst intervals with ΔtIBI < 1.0 s as not representative for single burst behavior. This cut removes “double-dipping” below our Γ threshold, caused for example by two subsequent activity “spikes” in the same burst with just enough recovery in between to get above threshold (very small ΔtIBI) or the interference of two burst waves (small ΔtIBI). We find the following grid averages: 〈ΔtIBI〉λ = (4.44 ± 0.26) s, 〈ΔtIBI〉λ = (4.0 ± 1.1) s, and 〈ΔtIBI〉λ = (3.5 ± 1.5) s. We see that the mean of ΔtIBI is decreasing roughly by 0.5 s per 3 mm reduction of λ; whereas the standard deviation increases considerably with decreasing λ, reflecting the more diverse spatial distribution of the burst patterns. According to these simulations, we can expect that spatial differentiation—the size of the burst suppression patterns and their timing—is intimately linked to the effective extent of the brain connectivity propagating the burst activity.

Figure 6

(A) Snapshot of the h activity of cortex at simulation time 9.54 s at 0.25 mM isoflurane with λ = λ1 = 2.7 cm. A black arrow shows the cortical location at which the corresponding time was recorded. An animation is provided as Movie 3 in the Supplementary Material (B) Snapshot at simulation time 5.79 s with λ = λ2 = 2.4 cm, the standard value. An animation is provided as Movie 4 in the Supplementary Material. The characteristic size of the burst patterns is reduced. (C) Snapshot at simulation time 11.70 s with λ = λ3 = 2.1 cm. An animation is provided as Movie 5 in the Supplementary Material. The characteristic size of the burst patterns is reduced even further. (D) Time series of h (blue) and Γ (green) taken from these three simulations, marked as “λ1,” “λ2,” and “λ3.” We see that inter-burst interval remains roughly the same.

4. Discussion

We find that the simulation of isoflurane induction with the model proposed here reproduces at least qualitatively the electrophysiological response that one can measure in the EEG (Foster et al., 2008). At light anesthesia there is an oscillatory shift to lower frequencies with higher amplitudes, then in deep anesthesia we find burst suppression patterns, and finally for even higher concentrations these bursts are abolished as cortex slowly heads toward electrotonic death. Significantly, in the burst suppression phase the prior regular activity of cortex is roughly “echoed” in the frequency content of the burst (Kroeger and Amzica, 2007; Ching et al., 2012), though of course the amplitude of the oscillations is much increased. We find that burst suppression results in dynamic and complex burst patterns that travel across cortex in waves, rather than remaining statically in place. The burst suppression phase is foreshadowed by the continuous elevation of peak PSP amplitudes, until finally these strong inputs destabilize the neural system into bursting. At maximum oscillation of the mean soma membrane of the neural population, strong depletion of the synaptic system leads to a sudden drop of the PSP peak amplitude, which suppresses the burst until the synaptic system is able to recover again. The relatively slow time scale of this recovery is what governs the periodicity of the bursts in the burst suppression regime in this model.

While we are mostly interested here in investigating the spatial differentiation of burst suppression qualitatively, the emergence of burst suppression in our simulations is also in rough quantitative agreement with what has been observed clinically. Because the emergence of burst suppression in the EEG represents a distinct endpoint, it has been proposed that it may be a suitable measure by which to titrate the administration of anesthesia to ensure optimal hypnosis. On this basis a variety of efforts have been made to estimate the concentration dependent emergence of burst suppression during anesthesia. It has been found that during sole agent isoflurane anesthesia, the burst suppression pattern can emerge at end-tidal concentrations as low as 1.2% (Hoffman and Edelman, 1995; Pilge et al., 2014), i.e., aqueous concentrations of ≃ 0.25 mM at 37°C. However when arterial blood concentrations of isoflurane have been measured the onset of burst suppression has been reported for levels as low as 34.9 μg/ml or ≃ 0.19 mM (Loomis et al., 1986) – close to the value at which we observed the onset of burst suppression in our model.

However, more importantly our model predicts the appearance of large-scale spatial burst patterns across cortex, which emerge, travel and disappear over time. In consequence, coherence of burst timing is mostly local, though one can expect to see characteristic burst onset time shifts in the case of burst patterns traveling across neighboring recording sites. Some of the simulated patterns become large in size intermittently, recruiting large parts of cortex and thus leading to more “global” correlations of burst timing. However, for the most part the complex spatiotemporal dynamics will lead to “local” correlations, with synchronization over large distances being mostly accidental. These predictions are at least qualitatively in line with the experimental observations of Lewis et al. (2013), see in particular their Figures 3, 4. We suggest that spatially dense experimental recordings may allow one to track such spatiotemporal burst patterns in detail. At least in principle it should be possible to reconstruct the underlying cortical state in terms of our model from such recordings, in particular if one tracks the activity over some length of time. The relatively slow spatiotemporal dynamics of the bursts may help in this regard.

We have shown here as well that increasing the inhibitory depletion factor f, or equivalently the resting values Γ and Γ of the inhibitory peak amplitudes, can abolish bursting in the model in a localized manner (in a patch of simulated cortical tissue). This is due to affecting the balance of two competing effects: the decrease of h reduces, whereas the increase of Γ increases, the capacity of cortex for self-excitation. For large enough increases of inhibition the former dominates the latter. This provides further possibilities for spatial differentiation in cortex as some tissue may have naturally higher inhibitory peak amplitude resting values, and hence be less capable of participating in burst suppression. In addition, it is expected that spatial heterogeneity in the cortical actions of anesthetics will contribute to the spatial differentiation of burst suppression. Most anesthetics that induce burst suppression are GABAergic agents, which in addition to enhancing inhibitory PSP action also produce increases in tonic inhibition, as well as reductions in tonic excitation by altering the activity of a variety of membrane bound channels that include two pore potassium channels, extrasynaptic GABAA and nicotinic acetylcholinergic ionotropic receptors. However, these synaptic and extra-synaptic channels, which exist in multiple isoforms that are variably affected by the same anesthetic agent, are not distributed uniformly throughout cortex. Thus, we would expect no two regions of cortex to share exactly the same propensity to burst for a given anesthetic level. This is not reflected in the current work, but left for future studies.

The hypothesized synaptic basis for spatially heterogeneous burst suppression suggests that if one can accordingly manipulate the cortical tissue, then one can artificially suppress, abolish or even enhance its participation in bursts. Importantly, our model predicts that drugs increasing inhibition could have the paradoxical effect of increasing burst activity, depending on the precise balance of the h decrease and Γ increase that they induce. In particular, one would typically expect that any paradoxical effects would occur at lower doses, since strongly increasing inhibition should eventually see the h decrease win over the Γ increase. It is interesting to note that a wide variety of GABAA modulators appear to have paradoxical effects at low doses, see for example (Bäckström et al., 2011) and references therein. Furthermore, our spatial model predicts that bursting activity of surrounding tissue can propagate into tissue that is incapable of bursting itself, to a depth depending on the density of synaptic connectivity, and lead there to “quasi-bursts” which simply reflect the dramatic variation of the synaptic input. However, if bursts are abolished by the mechanism suggested here, namely an increase in the inhibitory synaptic depletion factor f, then this tissue would show particularly large inhibitory PSP amplitudes and largely constant (rather than strongly varying) excitatory ones. This suggests that one could experimentally distinguish between bursting and “quasi-bursting” tissue by monitoring the size of the PSPs.

Furthermore, we have shown that both the spatial extent of the burst patterns, and the timing of the bursts (in particular the interval between bursts) depend on the characteristic scale of the brain connectivity effectively involved in propagating this activity: the shorter range these connections, the smaller the regions of coherent burst activity, and the more rapidly one burst follows on after the other. While at present this constitutes a qualitative finding, and while the implementation of brain connectivity in our model (homogeneous, isotropic and exponentially reducing with distance) is too simplistic to speak directly to the complexity of actual cortico-cortical connectivity, this nevertheless suggests that there is an intimate link between the spatiotemporal profile of burst suppression and the underlying brain connectivity. This will have to be taken into account when trying to improve the realism of such simulations. Thus, by being more specific about the anatomical structure of our mesoscopic model it may become possible that observations of burst suppression patterns will enable estimation of the effective connectivity of the bursting tissue. This would provide a new window on a difficult to access but key property of the brain.

Our present simulation has been restricted to studying the role that one slow modulatory system might have in the genesis of the burst suppression pattern. However, given the feedback inherent in the physiological and anatomical organization of cortex it is certainly only one of many systems that are capable of modulating cortical excitability, and hence the emergence of fast-slow bursting activity. Indeed we might hypothesize that such slow modulatory systems will span a number of functional scales in the brain—from perturbations in the autoregulatory coupling of neuronal and metabolic activity to alterations in the dynamics of cortico-thalamic and cortico-cortical feedback, to mention only the most obvious. However, regardless of the specifics of the slow system it is clear that any theory purporting to account for the genesis and features of cortical electrodynamics must be able to account for the reversible emergence of burst suppression in response to the action of anesthetic agents. In this respect both our model and the model of Ching et al. (2012) may be seen as meeting this requirement, even though they take as their starting points neuronal activity modeled at different spatial scales. A possible advantage of our approach, besides being able to deal easily with a spatially extended cortex, is that the modeled action of anesthesia is directly coupled to the emergence, and subsequent disappearance, of bursting. In contrast, in Ching et al. (2012) the parameters defining anesthetic action (τGABA and gGABA) are not directly related to the parameter JATP (the metabolic production rate of ATP) that defines the emergence of bursting.

Our model of synaptic depression is driven by pre-synaptic firing rates. This is unproblematic for the quasi-instantaneous local activity. However, an issue arises due to the propagation of action potentials from distant sources with finite velocities along cortico-cortical fibers. The current formulation of Equations (31, 33) as to how such distant inputs drive local PSPs is not entirely faithful to the actual physiology: the pre-synaptic firing rate is now modulated with the concurrent synaptic dynamics at the distant pre-synaptic site, which is then propagated with conduction delay; whereas it would be more physiologically accurate to propagate the pre-synaptic firing rates, and then modulate these conduction-delayed inputs with local synaptic dynamics. However, this would mean tracking separately local synaptic dynamics driven by quasi-instantaneous local (C) and delayed distant (CΦ) pre-synaptic firing rates, respectively. Our current model requires only one local synaptic dynamics (C), making it computationally simpler. Furthermore, one can argue that the current formulation better separates distant sources: If there are two distant sources, but only one of them begins to fire at higher rates, then only its signal would become depressed by the resulting synaptic dynamics, but not that of the other. Whereas if we were to drive a local CΦ directly with the sum of the delayed signals, then the higher firing of one would lead to synaptic depression also of the other. Yet one can argue to the contrary that such a conflation of distant inputs is simply part of the averaging approximation involved in considering neural populations rather than individual cells. We have made here a choice of convenience, but this may not always be possible. The parameters characterizing the (excitatory) synaptic dynamics are in this work taken to be spatially homogeneous. We are not aware of any extant empirical evidence that would require a significant deviation from homogeneity, but the alternative method mentioned above would be more suited for this case. For as currently formulated, the synaptic dynamics will be determined by parameters located at the origin of such activity and not, as must be physiologically the case, its pre-synaptic termination. A further possibility that we have not considered, due to the lack of any clear empirical data, is that these parameters themselves depend systematically on anesthetic concentration. If such a relationship is demonstrated, then this could act as a possible source of spatial inhomogeneity in the parameters of the synaptic dynamics, under the assumption that anaesthetic action shows spatial variability.

Finally, there are obvious extensions to this work that should be considered, but were beyond the scope of this initial investigation. The most obvious extension is to more systematically consider the effects of variations in the neural field model parameters defining the resting (unperturbed) EEG, and to study the resulting dynamics for spatially homogeneous models of synaptic resource depletion and anesthetic action. Furthermore, the neural field model used here can be extended to an equivalent neural mass mesh constrained by a real cortical head model based on MRI data (Bojak et al., 2010, 2011). Then one could investigate regional variations of the neural mass parameters with areal boundaries defined according to any of the available cortical structural/anatomical atlases [e.g., the Harvard-Oxford cortical and subcortical structural atlas (Desikan et al., 2006), or the Jülich histological atlas (Eickhoff et al., 2005)]. For example, each region could be assigned a distinct parameter set identified as producing physiologically plausible EEG within the physiologically admissible/plausible parameter space (Bojak and Liley, 2005). One could then additionally study regional variations in the synaptic resource depletion and/or anesthetic action parts of the model, as well as the interactions of heterogeneities in the component systems.

A perhaps less immediately obvious development, which however will be necessary to obtain a deeper understanding of the dynamical mechanisms responsible for the emergence of burst-suppression, would be some form of systematic bifurcation analysis. In such an analysis, the slow system would be “frozen”, i.e., one would set the synaptic ∂C/∂t ≡ 0, and a bifurcation analysis of the remaining “fast” subsystem would be performed by treating the C as bifurcation parameters. Such a bifurcation analysis, known as a fast-slow analysis, was pioneered by Rinzel (1985) in his formal analysis of bursting in biophysical models of the neuronal action potential. While such a bifurcation analysis is relatively straightforward for the temporal dynamics of non-linear ODE systems, using a variety of available bifurcation software tools like AUTO, Content, and MatCont (Meijer et al., 2013), it is considerably more challenging in systems of non-linear PDEs of the type we have studied here. A fast-slow bifurcation analysis of the spatiotemporal dynamics of burst suppression will require the development of new numerical methods and tools, which are only just beginning to emerge (Green and van Veen, 2014).

Our work then represents only a first step toward a deeper understanding of the spatiotemporal dynamics of burst suppression, in particular as induced by anesthesia. Yet it is already clear from the theoretical results obtained here, which were motivated by the recent experimental results of Lewis et al. (2013), that the classical understanding of burst suppression as spatially homogeneous phenomenon has become outdated. This can only add to the importance of burst suppression as a dynamical probe to investigate the properties and function of cortical tissue, whether in a theoretical modeling or an applied clinical setting. We expect that in the near future theory will be further challenged by the rapid technological advances in electrophysiology and neuroimaging, which are producing increasingly accurate and dense measurements of neuronal activity and cortical dynamics. This hopefully will allow us to test to what extent the synaptic resource depletion mechanism proposed here is indeed the driver of the observed burst suppression dynamics.

Author contributions

IB and DTJL developed the model. IB and ZVS programmed the simulation software, analyzed the generated data, and produced the figures and movies. IB and DTJL wrote the paper. All authors have approved the final version of the paper and the Supplementary Material.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.