General Anesthetic-induced Seizures Can Be Explained by a Mean-field Model of Cortical Dynamics

Many neurobiologic phenomena arise predominantly from interactions between and within populations of neurons, rather than directly from the properties of individual neurons. It is necessary to develop a theoretical framework in which to make sense of the experimental observations of these phenomena. What is the best model? Clearly there are many answers to this question. At one extreme are “bottom-up” models that concentrate entirely on biologic accuracy. These take the form of intensely detailed multicompartmental simulations of current flux in each channel in each neuron, and the exact neuronal morphology and connections. This may be termed the neuron-by-neuron approach. At the other end of the spectrum is the abstract neural network/computational approach. In these models, the only interest is in the flow of information, and they barely acknowledge any details of the real underlying biologic “hardware.” The mean-field model—a standard method of analysis in statistical physics—is a compromise between these two extremes. It assumes that the important aspects of neuronal function can be described by the average behavior of a small population of perhaps 10⁵neurons. Its advantages over the neuron-by-neuron methods are as follows:

1. Mean-field models are computationally manageable without a supercomputer cluster. In contrast, neuron-by-neuron methods for the cortex would involve solving approximately 10–100 sets of equations for up to 100 billion neurons, and a 100 trillion synapses in time steps of much less than 1 ms.

2. The simplification inherent in the mean-field model requires the modeler to try to extract the important components of the neuronal function and interactions on the basis of previous knowledge. This abstraction allows multiple redundant systems (which seem to be common in the central nervous system) to be simply modeled. For example, neuronal excitability may be increased by a large number of neuromodulators acting via  a plethora of mechanisms (e.g. , intracellular, ligand-gated, metabotropic receptors, or even via  the glia). In most cases, the end effect on the activity of the neuron is similar for all of these mechanisms and can be included in the mean-field model simply as a change in the sigmoid curve that describes the soma potential–to–firing rate relation—the average irritability of the neurons.

3. The output variable of the mean-field model is the mean soma potential of a population of neurons. This can be directly related to the easiest experimentally observed measure of cortical function—the electroencephalogram.

In the mean-field approach, the inputs and responses of individual neurons are replaced by population means. The mean-field model of the interactions between inhibitory and excitatory cortical neuronal populations has been previously used to describe the electroencephalographic effects of general anesthesia and sleep and has been described in detail in a recent article by Wilson et al.  7A synopsis is presented in the appendix to this article. It builds on a previous body of work by Liley, Wright, Robinson, Steyn-Ross, and others.8–13It has been used recently by Kramer et al.  14to study generalized seizure states in relation to epilepsy. We were at some pains to use realistic, physiologically derived values for the parameters of the model, as described in our previous articles. However, the question of accurate parametrization is still an open one because the measurements are typically derived from reported studies of single neurons, which are then assumed to be representative of average of the neural population. In brief, the model is formulated as a set of coupled differential equations in time and space that describe the time evolution of the mean soma potential of cortical “macrocolumns” of approximately 10⁵inhibitory and excitatory cortical neurons. The scalp electroencephalogram is an image of the time evolution of the summed mean soma potential of a large number of macrocolumns. In this way, the output from the theoretical model is a “pseudoelectroencephalogram” that can be directly compared with experimentally derived data that have been previously reported in the scientific literature. The differential equations also explicitly model action potential rates from cortical neuronal input. The distant input is modeled with a wave equation, and local input is modeled as an instantaneous feedback (no reverberating cortex–thalamus loops are included). The shapes and time durations of the excitatory postsynaptic potentials (EPSPs), and IPSPs, are modeled with a second-order differential equation in time. Reversal potentials are used to modulate the resulting input to the cell membrane, thus constraining the mean membrane potentials to lie within certain bounds. A sigmoidal function is used to relate mean action potential rate to mean membrane potential in both excitatory and inhibitory neuronal populations. Subcortical neuronal input is included through a random noise term. Specifically, we used these equations to model the cortex in two spatial dimensions as a 50 × 50-cm two-dimensional slab using periodic boundary conditions (approximately the area of unspecialized association cortex in humans). The quantitative effects of isoflurane and enflurane on IPSPs were incorporated using the data derived from Banks and Pearce4in a manner similar to that of Bojak and Liley.15 

Examples of the pseudoelectroencephalogram generated by the model in response to changes in enflurane concentration are shown in figure 1. In simple terms, the behavior of the model cortex shows three main patterns: (1) an active higher-firing state that is similar to that observed in the alert cortex (fig. 2, area A), (2) a quiescent low-firing state that is similar to that described in cortical neurons that are under the influence of common γ-aminobutyric acid–mediated anesthetic drugs (fig. 2, area B), and (3) a state in which the cortex is unstable and oscillates rapidly between silence and supramaximal firing. We identify this state with seizures (fig. 2, area C).

We can now vary the amplitude and time duration of the IPSP to mimic the effects of increasing concentrations of enflurane and isoflurane on our model cortex, and then see what changes in the state of the cortical firing occur. As shown in figure 2, the effects of the anesthetic drugs are visualized using a plot of the decay time constant of the IPSP (y-axis) versus  the time-integrated IPSP magnitude (x-axis). In essence, both enflurane and isoflurane cause an increase in total synaptic charge transfer (the area under each IPSP) of between 100 and 200%. This results in the cortex making a transition from an active state (fig. 2, area A)—which is identified as being the awake state—to a quiescent (anesthetized) state (fig. 2, area B).10,Figure 2bshows this transition more clearly by explicitly plotting the third dimension (the excitatory mean neuronal soma potential) along the z-axis. The trajectories followed by the model cortex in response to increasing concentrations of each anesthetic agent (enflurane = dashed line, isoflurane = solid line) are shown in figure 2. The start of each trajectory (bottom left) is zero anesthetic concentration, and the end of each trajectory (top right) is approximately 2 minimum alveolar concentration (MAC)⁵⁰. The point at which the anesthetic drug concentration equals 1 MAC⁵⁰is shown by X  s near the apex of each trajectory. It can be seen that both isoflurane and enflurane cause the cortex to move from the awake state (area A) to the anesthetized state (area B). However, there is a region (area C) where the cortex undergoes a sudden change in behavior to the unstable seizure state. Mathematically, the technical term for this is a subcritical Hopf bifurcation, to enter a limit-cycle oscillation. The figure shows that the enflurane trajectory has a greater probability of crossing, or entering, this unstable region, because enflurane reduces the IPSP magnitude more than isoflurane for a given IPSP time constant (i.e. , for an equal IPSP area, the shape of the enflurane-affected IPSP is “flatter” than that for isoflurane). When in this state (area C), the cortex alternates (approximately 5/s) abruptly between zero firing and supramaximal firing states, similar to that observed during seizures.16,17The changes in stability of the cortical model were initially derived from analysis of the eigenvalues of the equations and then checked by the computer simulations. In addition to states similar to spike-wave seizures, spiral waves were also sometimes observed in the model when the trajectory crossed an unstable region. The shape of the IPSP is critical. Therefore, for a relative IPSP charge-transfer of two, the cortex could be either in a state of coma, if the IPSP was short and of high amplitude (fig. 2, position shown by diamond and subgraph ii), or in a seizure state, if the IPSP was long and of low amplitude (fig. 2, position shown by square and subgraph i). The model also predicts the propensity of enflurane to induce seizures in almost all cases at higher concentrations of above approximately 2%—where the trajectory turns back into the middle of the unstable region—as the total IPSP charge decreases in the presence of a slow IPSP decay time. This is because the effect of enflurane to decrease the peak amplitude of the IPSP is greater than its effect in prolonging the decay time. These results are in agreement with experimental observation18and the theoretical work of Liley and Bojak.19 

Although we have chosen our model parameters of the figure to make the difference in stability of the two trajectories obvious, the parameters are all physiologically plausible, and the implication that enflurane has the greater propensity to cause seizures is clear. Also note the magnitude of the uncertainties in the Banks and Pearce data. To check how robust the model is to changes of parameters, we performed 1,000 runs with different, but plausible, parameter sets. We find that, for isoflurane, the mean percentage of points on the trajectories lying in an unstable region is 2.8–3.0%, with a total of 50–52% of trajectories crossing the unstable area at some point on their path. (The range gives an indication of the standard uncertainty in the result.) For enflurane, the same values are 11.5–12.7% and 57–60%. Moreover, 4.9–6.3% of trajectories for enflurane were unstable over at least half of their length, whereas no occasions of such a large instability were recorded for isoflurane.

Counterintuitively, the addition of thiopentone—usually an efficacious antiepileptic drug—to low concentrations of enflurane has been shown to increase seizure activity.20Thiopentone prolongs the IPSP decay time but has no effect on the peak amplitude of the IPSP21and therefore would generate a trajectory that travels up and rightward. From inspection of figure 2, it can be seen that the effect of adding thiopentone to a low concentration of enflurane could shift the enflurane trajectory up and to the right, thus causing it to be more likely to enter the oscillatory seizure state (the lower part of area C). Conversely, at higher concentrations of enflurane (> 2%), thiopentone has an anticonvulsant effect because it would shift the trajectory to the right, away from reentering area C.

The main interest of the model lies in its explanation of the proconvulsant effects of seemingly neurodepressant drugs. However, decreased excitation (modeled as shortening of the EPSP decay time constant and consequently a reduction in integrated magnitude) results in the area of instability (fig. 3A, area C) becoming narrower and moving upward and to the left. This explains the experimental observation that N -methyl-d-aspartate antagonist drugs effectively prevent or abort enflurane-induced convulsions.22Simulations of N -methyl-d-aspartate blockade with an EPSP reduced in both area and time by 25% give just 19–21% of 1,000 simulated enflurane paths as crossing the unstable region, compared with 57–60% of enflurane paths without the N -methyl-d-aspartate blockade.

Equivalently, the increased inhibition of neuronal firing caused by opening of extrasynaptic γ-aminobutyric acid/Cl channels or potassium channels may be modeled as hyperpolarization of the resting membrane potential. This also causes the size of area C of figure 3to decrease and shift to the left (fig. 3c). As an aside, it may be noted that both diminution of EPSP and resting membrane potential hyperpolarization also move the junction between the awake (area A) and anesthetized (area B) states to the left. This explains the fact that the anesthetic potency of volatile agents is increased by interactions with ketamine (an N -methyl-d-aspartate blocker) and with benzodiazepines (potent antiepileptic drugs that are known to increase tonic chloride currents).23 

Conversely, increased excitation induced by proconvulsant volatile “anesthetic"” agents such as 5-(2-cyclohexylideneethyl)-5-ethyl barbituric acid24(modeled as a 25% increase in the area of the EPSP) cause a great increase in the area of the unstable region (fig. 3b, area C), such that it dominates the phase plane. Likewise, decreased tonic inhibition—such as may be induced by physostigmine or 4-amino-pyridine25—may be modeled as depolarization of the resting membrane potential of the cortical neurons. This effect will greatly increase the propensity for the development of seizures (note large area C in fig. 3d).

It is easy to imagine that seizures could occur when the balance between inhibition and excitation in the cerebral cortex is unbalanced toward excessive excitation. It is less clear how seizures could occur when the cortex is in a predominantly inhibitory mode. Even this relatively simple model can quantitatively describe how it is possible for seizures to occur in an anesthetized cerebral cortex. It also shows that modest changes in the pattern of synaptic response can result in gross changes in the behavior of populations of neurons. The seizure state emerges as a natural consequence of the model dynamics. We contend that it is not possible to explain the phenomenon of seizures occurring in a quiescent cortex, without some attempt to model the interneuronal dynamical interactions. The critical factors are that (1) the system is in an inhibitory state and (2) there exists a disproportionate prolongation of the IPSP, which introduces a lag in the system, thus predisposing the cortex to oscillate. The impact of synaptic excitatory inputs on the generation of action potentials at the soma is controlled by complex active and passive, spatial and temporal filters in the dendritic tree.26,27It may be surmised that other drugs and pathologic processes that slow or broaden the dendritic summation will have similar proseizure effects.

The parameters that control our model (e.g. , the various voltages; the size, number, and time courses of the synaptic events; see appendix) are derived as far as possible from data that have been experimentally verified in the peer-reviewed literature. However, there exists some uncertainty as to the exact values. We justify our choice of parameter values from three arguments: (1) The experimental measurements of neuron-by-neuron parameters that are reported in the literature seem to vary widely with time and context. Therefore, narrow limits in the variability of individual parameters are probably less physiologically critical than relations between groups of parameters. (2) The fact that the cortex can continue to function, to some degree, in the face of significant neurobiologic disturbances (e.g. , histamine or muscarinic blockade by drugs) argues that multiple parallel redundant pathways exist. (3) The cortex exhibits some degree of adaptive self-organization, which would hugely constrain the effective size of the parameter space.

In our model, we do not specifically include subcortical interactions, except as a source of nonspecific activating input and via  the indirect effect of changes in neuromodulatory systems on the parameter values of the cortical neurons. There is strong experimental evidence that thalamocortical circuits are not essential for the development of seizures.17We also assume that the only effect of the anesthetic drugs is on the IPSP, and we model the cortex as a simple two-dimensional sheet, without taking into consideration three-dimensional layering or variations between different regions within the cortex. These and other issues related to seizure onset and termination (such as second-order changes in chloride and potassium balance across the neuronal cellular membrane) may be investigated in more detail in future models. For example, the exaggerated frontal electroencephalographic effects of common γ-aminobutyric acid–mediated anesthetic drugs would clearly indicate heterogeneity in parameters between frontal and occipital lobes. Similarly, most of the corticocortical axons interact with inhibitory interneurons and apical pyramidal cell dendrites in the superficial layers of the cortex (cortical layers 1 and 2), whereas axonal input from the thalamus occurs predominantly near the pyramidal cell bodies (cortical layers 3–5). However, even our simple model seems to have adequate explanatory power to describe most of the observed phenomena, because these phenomena (anesthesia and seizures) are global features of cortical function.

We conclude that the technique of mean-field modeling of the cerebral cortex is able to explain parsimoniously and elegantly a constellation of counterintuitive neuropharmacologic experimental data.

The model quantitatively describes the interactions between populations of excitatory pyramidal neurons (subscript e ) and inhibitory interneurons (subscript i ). We have used the convention of a→b, as meaning the direction of transmission in the synaptic connections is from the presynaptic nerve a to postsynaptic nerve b. The subscript sc  indicates random subcortical input that is independent of the cortical membrane potential. The variables V^a  and Φ^ab  vary in time and space. The time evolutions of the mean neuronal soma membrane potential (V^a ) in each population (excitatory [e ] and inhibitory [i ]) of neurons in response to synaptic input (ρ^a ψ^ab Φ^ab) are given by the following set of equations:

where τ^a  are the time constants of the neurons, ρ^a  are the strength of the postsynaptic potentials (proportional to the total charge transferred), ψ^ab  are the weighting functions that allow for the effects of reversal potentials and are described by the equation

and Φ^abare the synaptic input spike-rate densities (flux) that are described by the following equations:

where γ^ab  are the synaptic rate constants, N  ^αis the number of long-range connections, and N  ^βis the number of local, intramacrocolumn connections. The mean axonal velocity is given by ν, and the characteristic length (the length at which the connectivity decays to 1/e ) is 1/Λ ea . The spatial interactions among and within the macrocolumns are described by the two equations

The population firing rate of the neurons is related to their mean soma potential by sigmoidal functions

where θ^a  describes the inflection point voltage, and ς^a  describes the SD of the threshold potential. The parameters and ranges used in our simulations are shown in table 1(e  and i  refer to values assigned to excitatory and inhibitory cell populations).