Computational Aspects of Anesthetic Action in Simple Neural Models

Computational models of single neurons and networks of neurons were modified to permit modulation of the voltage-gated and ligand-gated ion channels by volatile and intravenous anesthetics. The effects of volatile anesthetics on neural activity were then evaluated as a function of anesthetic concentration and other model parameters. The impact on model behavior was assessed by computer simulation and, where possible, examination of the mathematical structure of the models using analytic techniques.

Four computational models of individual neurons and networks of neurons were examined. These were the Morris–Lecar model of oscillation in the barnacle muscle fiber, 49the Pinsky–Rinzel simplification 50,51of the 19-compartment model of the CA-3 hippocampal neuron by Traub et al. , 52the Golomb–Rinzel network model of the reticular nucleus of the thalamus, 53and the Wang–Buzsáki model of hippocampal interneurons. 54These models were selected because of their relative simplicity and to illustrate specific aspects of anesthetic action. The Morris–Lecar and Pinsky–Rinzel models were used to examine the behavior of individual neurons as voltage-gated channel activity is modulated by anesthetics. The Golomb–Rinzel and Wang–Buzsáki models were used to examine anesthetic effects on simple oscillating networks.

Each of the biophysical models presented is built around the formalism developed by Hodgkin and Huxley. 55Membrane potential V^mis obtained from the differential equation for membrane potential:where C is the membrane capacitance, I^appis the externally applied current, and I^ionis the sum of the currents resulting from the ion fluxes through the membrane ion channels. Essentially, this equation is a statement that the sum of the currents across the cell membrane must equal the externally applied current. Individual ionic currents whose sum comprises I^ionare modeled as a potential across a conductance. Thus, the individual currents Iⁱare expressed as where gⁱis the conductance of the channel and Vⁱis the potential due to the ion conducted by the i-th channel (reversal potential). The conductance is a highly nonlinear function of time and membrane potential. A general expression for gⁱis where ḡⁱis the maximal conductance of the i-th channel and is proportional to the channel density, pⁱand qⁱrepresent gates which must both be in a permissive state for current to flow, and lⁱand mⁱare the corresponding exponents which can be thought of as the number of gates in one ion channel. Both pⁱand qⁱrange between 0 and 1 to represent the fraction of gates in the permissive position. For both pⁱand qⁱthis fraction is typically governed by an equation of the form where α(V^m) is the voltage-dependent rate constant governing the rate that the gate transitions from the nonpermissive to the permissive state and β(V^m) is the voltage-dependent rate constant governing the rate that the gate transitions from the permissive to the nonpermissive state. Typically, α(V^m) and β(V^m) are highly nonlinear functions of membrane potential that must be empirically determined for each channel. For some systems calcium concentration can be a particularly important determinant of α and β. Temperature dependence of the model can be introduced through the system of equations defined by Eq. 4. This requires the use of a function φⁱ(T), such that where φⁱ(T) is often of the form where T is the new temperature, T⁰is the temperature for which the system was originally configured, and Q^10-iis the coefficient describing the increase in rate of the i-th process for a 10-degree increase in temperature. Thus, the rates, but not the steady state solutions of the system of equations defined by Eqs. 5 and 6will be altered by temperature.

The above system of differential equations is a general description of a single compartment that can be considered the fundamental building block for constructing biophysical neural models. Often, a single compartment is sufficient to capture the behavior of interest as in the original Hodgkin–Huxley model or the Morris–Lecar model described below. To capture more elaborate neural behavior, particularly when investigating an anatomically complex neuron with a heterogeneous distribution of ion channels, multiple compartments can be electrically coupled together as in the Pinsky–Rinzel model, described below, and the more complex models which motivated it. Networks of single or multiple-compartment neurons can be generated by coupling them with one or more excitatory or inhibitory synapses as with the Golomb–Rinzel and Wang–Buzsáki models, below. Typically, synaptic currents are modeled beginning with equations like those of Eqs. 2–4, in which the conductance is governed by the presynaptic potential and the kinetics of the specific neurotransmitter. The contribution of each neuron in the network to the synaptic current(s) must be considered, and the connection geometries governing neural interactions can range from all-to-all, to random, to highly specific.

The Morris–Lecar model ( Appendix A) 49is useful for demonstrating the effects on individual neural activity of anesthetic modulation of voltage-gated ion channels. Although rather simple, this model displays many behavioral modes seen in more complex neural models and, therefore, is often used to demonstrate these different modes of behavior. 56Moreover, its relative simplicity permits more analytical approaches. It consists of a single compartment containing a noninactivating high voltage activated (HVA) Ca²⁺channel, a K⁺channel, and a leak current.

The hippocampus plays an important role in memory and, because its anatomy is well suited for slice preparations, has been studied extensively in the laboratory. In addition, it has been the motivation for multiple computational modeling studies of increasing complexity, where this complexity extends to more realistic descriptions of ion channel activity, dendritic anatomy, distribution of ion channels over the dendrites, and network interactions. 57,58The Pinsky–Rinzel 50,51( Appendix B) reduction of Traub et al.’  s 19-compartment hippocampal neuron model 52is used to illustrate the consequences of volatile anesthetic modulation of ion channel activity in hippocampal neurons. The two-compartment Pinsky–Rinzel model consists of single somatic and dendritic compartments that are electrically coupled. The somatic compartment contains a fast Na⁺channel, a fast K⁺channel, and a leak current. The dendritic compartment contains an HVA Ca²⁺channel, two Ca²⁺activated potassium channels (BK and SK), and a leak current.

Although the thalamus is primarily known for its multiple nuclei relaying afferent sensory input, the reticular nucleus of the thalamus is thought to be the source of the spindle oscillations (7–14 Hz) seen during drowsiness, at sleep onset, and barbiturate anesthesia, and is the source of the electroencephalographic α rhythm (8–12 Hz). 43,47,48The behavior of this nucleus is thought to be modulated by excitatory input from thalamocortical neurons which, themselves, receive several types of inhibitory input from the neurons of the reticular nucleus. Computational models of this system vary in complexity with respect to the description of the ion channels of the individual neurons as well as the nature and complexity of the interactions of the reticular nucleus with thalamocortical neurons. 53,59,60A simplified model of the reticular nucleus of the thalamus was developed by Golomb and Rinzel ( Appendix C) 53in which membrane potential is controlled by a low voltage activated (LVA) Ca²⁺channel and a leak current, and neural interactions are mediated by GABA^Ainhibitory interactions. Since modulation of GABA^Areceptor activity is thought to be the basis for the action of a number of anesthetics, this model provides the opportunity to computationally examine the effects of anesthetics on a physiologically important neural system that may be related to consciousness. The model in  Appendix Chas been modified slightly from the original to permit partial connectivity between neurons in the network.

As delineated in the Introduction, interactions of the neocortex may be regulated by a network of fast-spiking interneurons. These fast-spiking interneurons may also contribute to the behavior of the hippocampus. 61,62Wang and Buzsáki (see  Appendix D) 54have constructed a model of this network in which neural interactions are modulated by GABA^Ainhibitory interactions. Because the voltage-gated ion channels that define the properties of the individual neurons in this network are different from those of the reticular nucleus of the thalamus, the network properties also differ. This difference arises primarily because the output of the thalamic neurons described above is a burst controlled by Ca²⁺entry as opposed to spikes controlled by the dynamics of fast-sodium and fast-potassium channels.

The volatile anesthetics are known to modulate both voltage-gated 63and ligand-gated 64ion channels. Of the voltage-gated ion channels, the effect on the Ca²⁺channels is most pronounced, 63although other types of receptors may play a role. 65–69Inhibitory effects on both the LVA and HVA Ca²⁺channels decrease Ca²⁺influx under general anesthesia. 63Effects on the T- (LVA), L- (HVA), and N- (HVA) type Ca²⁺channels have been demonstrated. 63,70–72The decrease in Ca²⁺current as volatile anesthetic concentration is increased appears to be relatively voltage insensitive. 63Consequently, as a first approximation which neglects transient effects upon activation of the channel, 71it is reasonable to model volatile anesthetic modulation of Ca²⁺channels as a decrease in the maximal Ca²⁺conductance (ḡ^Cain Eqs. A.1, B.2b, and C.1) which varies with volatile anesthetic concentration. The extent that a channel is inhibited as anesthetic concentration is increased is frequently modeled with the Hill equation where [A] is the anesthetic concentration, n is the Hill coefficient, and EC⁵⁰is the anesthetic concentration at which channel activity is inhibited by 50%. This leads to the following modification of Eq. 3:where, effectively, the maximal conductance is now modulated as a function of anesthetic concentration. Volatile anesthetic effects on ion channels other than the Ca²⁺channels do occur. However, the effects on the fast-sodium and fast-potassium channels occur at much higher concentrations than for the Ca²⁺channels. 63,65,66Meaningful volatile anesthetic effects on the Ca²⁺-activated potassium channels may occur at anesthetic concentrations close to the EC⁵⁰for the Ca²⁺channels. 73–76 

A variety of intravenous and volatile anesthetic effects are known for the ligand-gated ion channels. Most of these involve the GABA^Ainhibitory channel that controls neural interactions in the Golomb–Rinzel model 53of the reticular nucleus of the thalamus and the Wang–Buzsáki 54model of the fast-spiking interneuronal network. The GABA^Ainhibitory channel is the only ligand-gated ion channel in these models.

The intravenous anesthetics such as the barbiturates, 77,78etomidate, 79propofol, 78,80and neurosteroids 81are all known to increase the inhibitory chloride current by prolonging the open time of the GABA^Achannel when it is stimulated with GABA. 80,82At higher concentrations of the intravenous anesthetics a decrease in the chloride current through the GABA^Achannel is observed for at least several classes of the intravenous anesthetics which are known to prolong GABA^Aopening. 83,84At somewhat higher concentrations the intravenous anesthetics listed above can activate chloride currents in the absence of GABA. 80,81,85,86The benzodiazepines, though exerting their action at the GABA^Areceptor, differ from the intravenous anesthetics in that they appear to exert their effect by increasing the rate at which the GABA^Achannel opens, 87and are also distinguished by their inability to directly induce chloride currents through the GABA^Areceptor. The effects of volatile anesthetics on the GABA^Achannel share several similarities with the intravenous anesthetics. 64Their effect is characterized by prolongation of channel open time. However, at somewhat higher volatile anesthetic concentrations, the amplitude of the chloride current in the GABA^Achannel is also decreased.

The above description of anesthetic effects on the GABA^Achannel leads directly to a set of modifications to the descriptions of GABA^Adynamics in the above computational models. An increase in the open time of the GABA^Achannel, as seen with the intravenous and volatile anesthetics, can be modeled by decreasing the rate constant which governs the rate of channel closure in the equations describing GABA^Achannel dynamics (β^synin Eqs. C.3 and D.4of the Appendices). This decrease in the rate constant as a function of anesthetic concentration has been measured experimentally for propofol and thiopental, 78and isoflurane, enflurane, and halothane. 64An increase in the rate of GABA^Achannel opening, as seen with the benzodiazepines, can be modeled by increasing the parameter that governs the rate at which the channel opens (α^synin Eqs. C.3 and D.4of the Appendices). A decrease in the amplitude of the inhibitory chloride current generated when the GABA^Ais open can be modeled by decreasing the synaptic conductance (ḡ^synin Eqs. C.1 and D.1of the Appendices). For some anesthetics, 64,78detailed quantitative descriptions of how these parameters vary with anesthetic concentration have been obtained and can be incorporated in the mathematical description of the model.

The models presented in the Appendices were implemented with Mathematica 4.02.2 88and were solved with the Runge–Kutta 89integration option. Additional simulations were also performed with the neural simulation language Genesis 902.0 that uses the exponential Euler integration scheme, where step size was progressively decreased and the stability of solutions with respect to chosen step size confirmed. All simulations were carried out until steady state behavior was achieved, and results from the terminal portion of the simulations are usually presented.

Analytic approaches become less practical with increasing model complexity. Consequently, these are considered only for the Morris–Lecar model. The example given here follows a very typical approach for problems of this type, 6and was previously exemplified for the Morris–Lecar model. 56Here, this approach is adjusted to examine the alteration in model behavior as anesthetic concentration is varied.

Although more complex questions can be addressed, the most fundamental issue is when and how qualitative behavior changes as one or more parameters are varied. In this case, the parameter is, effectively, anesthetic concentration through its modulation of ḡ^Ca(Eq. A.1in the Appendices). The first step is to find the equilibrium point(s) of the system by setting the derivatives of Eqs. 1 and 4to zero and solving for the variables V^mand pⁱ. For each of the equilibrium points, the system is linearized and the stability of the linearized system examined. This is typically done by computing the eigenvalues of the matrix which characterizes the linearized system at the given equilibrium point. If the eigenvalues of the matrix are given as λⁱ, this corresponds to modes xⁱ(t) of system behavior in the vicinity of the equilibrium point of the form

which indicates that the solution will tend to leave the vicinity of the equilibrium point for Re(λⁱ) > 0, will head toward the equilibrium point for Re(λⁱ) < 0, and exhibit a transition (bifurcation) in behavior for Re(λⁱ) = 0, where Re(·) denotes the real portion of the possibly complex eigenvalue. Thus, by examining λⁱas parameters of the system are varied important information can be obtained about the ability of given parameters to impact on fundamental aspects of system behavior. A typical behavioral change indicated by Re(λⁱ) passing through zero would be a switch from quiescence (Re(λⁱ) < 0) to sustained oscillatory behavior (Re(λⁱ) > 0).

All analytic calculations were carried out using both the symbolic and numeric capabilities of Mathematica. The equilibrium points were obtained numerically, these results were substituted into the corresponding Jacobian matrix, obtained symbolically from the linearized system, and the eigenvalues of this matrix were then computed.

When parameterized as in  Appendix A, the Morris–Lecar model produces sustained oscillations for injection currents from 33 to 42 μA/cm². In figure 1, for a fixed injection current of 35 μA/cm², the maximal Ca²⁺conductance (ḡ^Cain Eq. A.1) is decreased, as might occur in the presence of a volatile anesthetic, until oscillations cease. Initially, as ḡ^Cais decreased, there is a graded decrease in the frequency of the oscillations (not shown). As illustrated in figure 1, for the given injection current, only a 10.8% decrease in the Ca²⁺conductance is sufficient to prevent sustained oscillations, and the termination of oscillatory behavior occurs abruptly as ḡ^Cais decreased. The loci of injection current and ḡ^Cawhere oscillations cease is shown in figure 2A. Note that when ḡ^Cais decreased from its baseline value, larger injection currents are required for sustained oscillations. Conversely, larger injections currents will hyperpolarize the neuron unless ḡ^Cais decreased. The Morris–Lecar model was developed for a physiologic preparation maintained at 22°C (room temperature). Volatile anesthetic sensitivity of HVA calcium channels at room temperature is known, and demonstrates a midpoint of approximately 0.85 mm and a Hill coefficient of 1.5 for halothane. 63Using these values to parameterize Eq. 8, the data of figure 2Acan be used to construct the loci of volatile anesthetic concentration and injection current, shown in figure 2B, for which oscillations cease. Simultaneously incorporating volatile anesthetic sensitivity of the potassium channels by assuming that halothane suppression of potassium channel activity has a midpoint of between 3.0 63and 6.4 mm 66and a Hill coefficient of 1.5 10introduces no additional behavior. For the points shown in figure 2, the termination of oscillatory behavior as ḡ^Cais decreased corresponds analytically to a change in sign of the real component of one of the corresponding eigenvalues from positive to negative, demonstrating a basis for the abruptness of the transition.

The response of the Pinsky–Rinzel model of the CA3 hippocampal neuron (see  Appendix B) to steady current injection to the somal compartment is known to vary qualitatively with the magnitude of the injection current, 50as summarized in figure 3. Low levels of injection current lead to a bursting pattern (fig. 3A) mediated by calcium entry. Higher currents lead to complex patterns of spiking and bursting due to interactions between the somal and dendritic compartments (fig. 3B). Even higher currents lead to a somatic spiking pattern (fig. 3C), whose rate increases with increasing somal current injection (not shown) until the current is so great that the neuron remains in the depolarized state and cannot generate spike activity.

Halothane modulation of maximal Ca²⁺conductance (ḡ^Cain Eq. B.2b) was modeled using the same scheme and parameters as above. Modulation of the bursting pattern of figure 3Aby halothane is shown in figure 4for halothane concentrations of 0.10, 0.20, and 0.33 mm. Although a small increase in burst duration and frequency is seen in figure 4Afor concentrations of 0.10 mm, most of the increase occurs between 0.15 mm and 0.2 mm (fig. 4B). Between 0.2 mm and 0.33 mm, irregular behavior composed of shorter bursts and spikes eventually gives way to the regular spiking behavior shown in figure 4C. For somal injection currents that produce the irregular pattern of figure 3B, halothane modulation of Ca²⁺conductance produces a spiking pattern as halothane concentration is increased (not shown). Once the somal injection current is high enough to produce pure spiking behavior as in figure 3C, halothane modulation of Ca²⁺-conductance has no impact on the pattern and frequency of the spike train.

The alterations in burst duration and frequency occur as a direct result of decreases in the Ca²⁺current. Although a decrease in Ca²⁺current requires a smaller dendritic K⁺current to repolarize the neuron, the burst duration still increases for a range of volatile anesthetic concentrations because decreases in Ca²⁺accumulation delay production of sufficient dendritic Ca²⁺-dependent K⁺currents (primarily I^K-Cof Eq. B.2a) to terminate the burst. The decreased rate of intracellular Ca²⁺accumulation is due to both the decrease in the amplitude of the Ca²⁺current and, indirectly, to the corresponding decrease in the membrane potential of the dendritic compartment, which then prevents full activation of the gating parameter of the Ca²⁺channel (s in Eqs. B.2b, B.7 and B.8). The decrease in Ca²⁺concentration achieved during a burst also contributes to the observed increase in burst frequency. A burst cannot be initiated until the dendritic K⁺current (primarily I^K-AHPof Eq. B.2a) is sufficiently small, which occurs when the gating variable q for I^K-AHPdecreases below a specific threshold. When the Ca²⁺current is decreased by a volatile anesthetic, the decreased Ca²⁺concentration decreases the activation variable α^q(Eq. B.14a) and, hence, the maximal value of q achieved during a burst. Since the decay of q is relatively invariant with respect to membrane potential and Ca²⁺concentration (Eq. B.14b), a smaller maximal value of q will require a shorter interval to reach the threshold where burst initiation becomes possible, leading to an increase in burst frequency. At higher volatile anesthetic concentrations, the Ca²⁺conductance is diminished to the point where Ca²⁺mediated voltage spikes do not occur in the dendritic compartment. Burst activity becomes impossible and is supplanted by rapid somatic spiking, effectively switching the behavioral mode of the neuron.

Incorporating volatile anesthetic effects on the sodium and potassium channels in the somal compartment adds little new behavior. The suppression of sodium channel activity by halothane has been reported to have a midpoint that ranges from 2.0 66to 2.6 63mm, and that for the potassium channel from 3.0 63to 6.4 66mm. Following earlier work, a Hill coefficient of 1.5 is assumed. 10When the more sensitive value of 2.0 mm is used for the sodium channel, the anesthetic concentration at which regular spiking behavior emerges will increase by about 0.01 mm, but only if the EC⁵⁰of the potassium channel is close to 6.4 mm. Once regular spiking appears, it persists with a decreasing amplitude as anesthetic concentration is increased to at least several times the value of the midpoint of the sodium channel sensitivity.

The behavior of the Golomb–Rinzel model of the reticular nucleus of the thalamus as inhibitory interactions are modulated (fig. 5) for a network where, on average, each neuron interacts with 80% of the other network neurons. Here, the constant governing the rate at which the GABA^Achannels close is varied. This might occur in the presence of intravenous and volatile anesthetics that modulate GABA^Achannel behavior by prolonging the time that the GABA^Achannel remains open. Smaller rate constants correspond to more prolonged channel opening. Measurements of the rate constant β^synyield a value of 0.067 · ms⁻¹, which decreases to about 0.014 · ms⁻¹in the presence of 200 μm thiopental and 10 μm propofol. 78The volatile anesthetics decrease β^synfrom about 2.5 (halothane) to 5 (isoflurane) times its baseline value. 64,91As β^synis decreased and GABA^Achannels remain open longer, both the frequency of the overall neural activity and the coherence of this activity are affected (figs. 5, 6). The increased coherence as β^synis decreased is reflected in the increased amplitude and decreased frequency of the average membrane potential, which should be interpreted as a field potential analogous to the electroencephalograph. In contrast to overall network behavior, the gross behavior of individual neurons is minimally affected by changes in β^syn, as seen in figure 6.

Decreases in maximal Ca²⁺conductance (ḡ^Caof Eq. C.1) of about 20–40%, as might occur in the presence of the volatile 63or intravenous anesthetics, 92permits network synchrony to occur at somewhat larger values of β^syn(not shown). The network behavior that is most prominent with decreased Ca²⁺conductance is the increase in the value of β^synwhere the fundamental frequency of the average membrane potential approaches that of the individual neurons, signifying relatively complete network synchrony. This occurs for β^syn=.01 · ms⁻¹when ḡ^Ca= .4 mS/cm², and β^syn= .02 · ms⁻¹when ḡ^Ca= .3 mS/cm². The decreased Ca²⁺current permits a relatively greater impact of the inhibitory currents on membrane potential, effectively increasing the coupling between neurons, permitting synchronization to occur at larger values of β^syn. Unlike the Pinsky–Rinzel model of the CA3 hippocampal neuron presented above, the relatively simple calcium dynamics of the Golomb–Rinzel model of the reticular nucleus of the thalamus precludes refinement of the calcium burst by delayed calcium entry. Consequently, the shape of the calcium burst and the burst frequency remain unaffected with the above decrements in Ca²⁺conductance.

Both the volatile 64,91and intravenous 83,84anesthetics may decrease the amplitude of the inhibitory currents produced by GABA^Achannels. As might be anticipated, 53the above results are affected minimally by decreases of 10–20% in the conductance governing GABA^Asynaptic interactions (ḡ^synin Eq. C.1). Decreases of this magnitude might be observed at relatively high concentrations of the volatile anesthetics. 64An increase in the rate constant governing GABA^Aopening (α^synin Eq. C.3), as is hypothesized to occur in the presences of benzodiazepines, 87has little effect on the behavior of individual neurons or the entire network. This occurs because the onset of inhibitory currents is already very rapid when compared to the relatively slow time course of the individual neurons. Consequently, increasing α^syndoes little to alter the time course or intensity of inhibitory interactions.

As emphasized by Wang and Buzsáki 54and illustrated in figure 7, the synchrony of the fast-spiking interneuronal network depends critically on the rate constant governing the closure of the GABA^Asynapses. However, in marked contrast to the network model of the reticular nucleus of the thalamus, reduction in the value of β^synleads to decreased  synchrony of the network. This is reflected by an increase in frequency and decrease in amplitude of the membrane potential when it is averaged over all neurons in the network as a network analog of the electroencephalograph. Importantly, this occurs for relatively modest reductions in β^syn, especially when compared to the reductions in β^synthat were necessary to achieve network synchrony in the model of the reticular nucleus of the thalamus. The reduction in synchrony with reduction of β^synwas seen for a range of values of ḡ^syn. In contrast to the properties of the Golomb–Rinzel model depicted in figure 6, the frequencies at which individual neurons in the network oscillate decrease somewhat with decreasing β^syn.

The sensitivity of the fast-spiking interneuronal network to alterations in the rate at which the GABA^Achannels open, as might occur in the presence of the benzodiazepines, is shown in figure 8. An increase in the rate at which GABA^Achannels open can lead to greater degrees of synchronous behavior. Although the values of α^synused to generate figure 8are somewhat smaller than that of 12 · ms⁻¹, used to generate the simulations of figure 7, they demonstrate a sensitivity of the fast-spiking interneuron model to changes in the rate of GABA^Achannel opening that are not present in the more slowly oscillating model of the reticular nucleus of the thalamus.

Incorporating volatile anesthetic effects on Na⁺and K⁺conductance as described for the model of the hippocampal neuron has minimal impact on the behavior of individual neurons within the network. However, the decreased spike amplitude as anesthetic concentration approaches the EC⁵⁰of the Na⁺channel leads to decreased inhibitory interactions and loss of network synchrony for halothane concentrations of about 1.5 mm.

To help conceptualize the integrated response to anesthetic modulation of voltage-gated and ligand-gated ion channels and the contribution of this integrated response to general anesthesia, four computational neural models were modified to incorporate known features of anesthetic action on ion channel behavior. Models were chosen for their parsimony and their ability to illustrate integrative aspects of anesthetic action in systems with features that may be relevant to general anesthesia.

The Morris–Lecar model of the barnacle giant muscle fiber, 49although of minimal clinical relevance, is an important model of calcium oscillations whose simplicity permits some degree of mathematical analysis which complements observations made from the simulations. The major finding with this model is that relatively modest decrements in maximal Ca²⁺conductance, as might occur in the presence of the volatile anesthetics, can terminate the oscillations. Importantly, the Ca²⁺conductance and, therefore, the anesthetic concentration at which neural quiescence emerges is a function of the level of stimulation provided by the amplitude of the injection current, consistent with the observation that general anesthesia is a function of the level of stimulation. An examination of the underlying mathematical structure of the model reveals not just a change in the nature of the oscillation, but the emergence of a perfectly stable system that is no longer capable of oscillation.

Although it is known that volatile anesthetics can suppress Ca²⁺current, obtaining model behavior as a function of the concentration of a particular volatile anesthetic, as was done in figure 2for the Morris–Lecar model, is still somewhat imprecise. First, the kinetics of the volatile anesthetic's interaction with the Ca²⁺channel are not well-established, nor is the effect of the volatile anesthetic on the kinetics of the Ca²⁺current. Second, the full complement of necessary experimental data are unavailable. Minimally, this would include measurements of anesthetic modulation of ion channel activity for the same type of preparation as the model, obtained at the same temperature as the experimental data on which the model is based, and with sufficient measurements to describe channel behavior over the full range of anesthetic concentrations. The Morris–Lecar model was developed for a preparation of the barnacle giant muscle fiber at room temperature (22°C). Data obtained at room temperature from a different type of cell demonstrates values of 0.85 mm for the EC⁵⁰, and 1.5 for the Hill coefficient (n in Eq. 7), 63which were used along with Eq. 7to estimate volatile anesthetic effects for other anesthetic concentrations. Even though there is some imprecision in relating volatile anesthetic concentration to decreases in calcium current, this does not alter the fundamental observations of threshold-like behavior occurring for decreases in maximal Ca²⁺conductance well below 50%, and that the threshold varies with the level of stimulation current.

The Pinsky–Rinzel 50model of the CA3 hippocampal neuron was examined with the same computational approach used to examine the Morris–Lecar model of the barnacle giant muscle fiber. To determine the impact of volatile anesthetics on model behavior, the model incorporated the same quantitative relationship between anesthetic concentration and maximal Ca²⁺conductance (ḡ^Cain Eqs. B.2b and B.8) used above in the Morris–Lecar model. Volatile anesthetic modulation of Ca²⁺current leads to an increase in burst frequency, a lengthening of burst duration and, finally, a switch from bursting to spiking behavior as anesthetic concentration is increased (fig. 4). As detailed in the Results, all of these effects occur in the model as a direct result of decreases in the Ca²⁺current and the subsequent decrease in intracellular calcium concentration. The increase in burst duration is a consequence of the decrease in the calcium-activated potassium current (I^K−Cin Eq. B.2a, also known as BK). Therefore, any decrease in this current by direct action of volatile anesthetics, 73–76should lead to additional burst duration. Because the increase in burst frequency is tied to reductions in the afterhyperpolarization current (I^K-AHPin Eq. B.2a, also known as SK), it can be seen that if volatile anesthetics had the additional effect of reducing the afterhyperpolarization current, then burst frequency might be further increased.

Again, these results demonstrate that relatively modest volatile anesthetic effects on voltage-gated ion channels could have important consequences for the behavior of the neuron. The prolongation of burst duration in the presence of volatile anesthetics is known to occur in vitro . 93A switch from a bursting mode to a tonic-spiking mode for isolated hippocampal neurons in the presence of volatile anesthetics is not as well established experimentally, 94where network considerations may predominate. However, in simpler systems, a sufficiently high Ca²⁺conductance is necessary for bursting behavior to supplant spiking behavior. 95The ability of a volatile anesthetic to switch a neuron from a bursting pattern to a tonic spiking one whose rate is a function of the input current could have important consequences for information transfer in the brain. Switching of neural firing patterns is already used by the brain to regulate information flow through the thalamus where modulation of a T-type Ca²⁺channel determines the mode of behavior, 96and bursting is more likely to elicit cortical activity. 97 

As with the Morris–Lecar model, the relationship between specific effects in the model and a given anesthetic concentration is imprecise for many of the same reasons. An additional issue with the model of the CA3 hippocampal neuron used here is the presence of only a single type of Ca²⁺channel, when multiple HVA and LVA Ca²⁺channels contribute to the membrane potential of CA3 neurons, 58each with a unique sensitivity to the volatile anesthetics. 63,70–72Ultimately, the full impact of the effects of volatile anesthetics on the model neurons may not be realized until they are incorporated in a network which also incorporates the volatile anesthetic effects on inhibitory interactions 64as well as the contributions to network dynamics from other inhibitory and excitatory interactions which may not be affected by the volatile anesthetics.

An appreciation for the contribution of network interactions can be seen in the simulation of the Golomb–Rinzel model of the reticular nucleus of the thalamus. Here, a decrease in the rate at which the GABA^Achannels close, which leads to prolongation of the inhibitory current, was used to incorporate the effects of intravenous or volatile anesthetic effects on inhibitory synaptic interaction. The anesthetic-induced reductions in the rate of GABA^Achannel closure (β^synin Eq. C.3) that are introduced in the model correspond to those found experimentally for clinically relevant concentrations of anesthetics. 64,78As the open time of the GABA^Achannels increased, progressive slowing of overall network activity, but not that of individual neurons, was observed, reflecting a growing synchrony in network behavior. Although there is a level of conductance for the GABA^Achannel (ḡ^synin Eq. C.1) below which synchrony becomes impossible, 53for the most part decreases in this conductance only altered the range of the rate parameter governing closure of the GABA^Achannels (β^synin Eq. C.3) over which the transition to synchronous behavior is observed. As emphasized above for the model of the CA3 hippocampal neuron, volatile anesthetics may also affect the currents of the voltage-gated ion channels, particularly Ca²⁺channels. Although the model of the reticular nucleus used here incorporates a T-type Ca²⁺channel, many model features that might permit more elaborate interactions of volatile anesthetics with the individual neurons, particularly the Ca²⁺-dependent potassium currents were not present. Despite this, decreases in Ca²⁺current increased the ability to achieve network synchrony because this permitted a relatively greater influence of the inhibitory synaptic currents on membrane potential. Increasing the rate of GABA^Achannel opening, as might occur with the benzodiazepines, did not affect model behavior since the time course of the opening was already quite fast relative to the time course of the membrane potential.

Additional aspects of anesthetic action on the behavior of the reticular nucleus of the thalamus could emerge with a more realistic model. Excitatory input from thalamocortical neurons, which are themselves inhibited by the reticular nucleus, is thought to contribute to regulation of the rhythm of the reticular nucleus. 48In addition, both  GABA^Aand GABA^Binhibitory interactions are present within the reticular nucleus and the inhibitory pathway from the reticular nucleus to the thalamocortical network. Moreover, the oscillations of both the thalamocortical neurons 98and the neurons of the reticular nucleus 99are driven in part by T-type Ca²⁺-currents which are known to be affected by volatile anesthetics, 70and whose modulation increased the synchrony of the model of the reticular nucleus, as shown above. Although GABA^Breceptors are not affected by volatile or intravenous anesthetics, 100the longer time constants associated with their activation and inactivation 101could impact on model behavior. When thalamocortical neurons and GABA^Binhibitory interactions are incorporated in models of the reticular nucleus, 59,60,102they appear to confer additional stability to the resulting oscillations. 59 

The Wang–Buzsáki model of a fast-spiking interneuron network, when stimulated with a level of input current sufficient to produce oscillations in the γ range (40–80 Hz), is notable for its paradoxical behavior when contrasted with the behavior of the reticular nucleus of the thalamus. Decreases in the rate constant governing closure of the GABA^Achannels (β^synin Eq. D.4), as occurs with many intravenous and volatile anesthetics, lead to a loss  of synchrony (fig. 7), where the theoretical 103and computational 54basis for this behavior have been further elucidated elsewhere. Increases in the rate constant governing opening of GABA^Achannels (α^synin Eq. D.4), which occurs with the benzodiazepines, may be more important to the behavior of this network (fig. 8) than in the relatively slower thalamic network. A much more complex model of γ oscillations in a hippocampal slice preparation, 8which incorporates the interaction of detailed multicompartment models of hippocampal neurons modulated by fast-spiking interneurons, requires more than a decrease in the rate of GABA^Achannel closing to lose the synchrony of the γ frequency oscillations. An additional GABA^Ainhibitory leak current, as is documented to occur at higher concentrations of the intravenous anesthetics, 80,81,85,86must be incorporated to induce the loss of synchrony of the γ frequency oscillations. The addition of gap junctions to permit electrical coupling between the dendrites of interneurons is a model enhancement that can affect synchrony in certain types of networks. 104 

The specific molecular site(s) responsible for general anesthesia remain(s) controversial. Since many of the intravenous and volatile anesthetics are known to exert their effect at the GABA^Areceptor, it is tempting to implicate these actions at the GABA^Areceptor as the unitary mechanism of general anesthesia. However, the volatile and intravenous anesthetics also affect voltage-gated ion channel activity, though it has been postulated that the anesthetic concentration at which these effects occur and the nature of these effects render them clinically irrelevant. 10The results presented above from the simulations with the Morris–Lecar model of the barnacle giant muscle fiber and the Pinsky–Rinzel model of the CA3 hippocampal neuron call this view into question. In each case, relatively modest reduction of Ca²⁺current corresponding to clinically relevant anesthetic concentrations brought about dramatic alterations in model behavior. Nevertheless, fundamental changes in network synchrony were observed when clinically relevant concentrations of intravenous or volatile anesthetics decreased the rate of GABA^Achannel closure, and network synchrony could be further modulated by concurrent effects at the voltage gated channels.

The behavior of all four models is summarized in figure 9, where the concentration-effect curves for halothane on the voltage- and ligand-gated ion channels are shown along with the halothane concentrations at which model behavior is fundamentally altered. Although relatively high concentrations of halothane are necessary to achieve network synchrony in the thalamic model, this may be consistent with known differences of halothane on the electroencephalograph compared with other volatile anesthetics. For example, if it is assumed that the maximal reduction in the rate of GABA^Achannel closing is twice as great for isoflurane as for halothane, 64,91then network synchrony is achieved for isoflurane concentrations of 0.40 mm (as opposed to 0.70 mm) in the thalamic model, and synchrony is lost for isoflurane concentrations of 0.16 mm (as opposed to 0.26 mm) in the fast spiking interneuron network.

General anesthesia is not a single dichotomous state, and includes sedation, amnesia, loss of consciousness, immobility, and blunting of autonomic reflexes, with each of these features appearing once anesthetic concentrations reach a given threshold. 18–22Threshold behavior appears in the models presented here as well as in a more abstract model of general anesthetic action, 7and was the criterion used to generate figure 9. Clinically and encephalographically, 25general anesthesia can be biphasic, with an initial excitement phase followed by increasing depth of anesthesia as anesthetic concentration is increased. Several of the models presented here exhibit behavior consistent with inhibition and several models exhibit behavior consistent with excitement. In the CA3 hippocampal model, the increase in burst frequency, increase in burst duration and, perhaps, the shift to a tonically spiking pattern of behavior as anesthetic concentration is increased could be considered examples of increased excitability. The loss of synchrony in the fast-spiking interneuron model as GABA^Achannels remain open longer could lead to less organized and, therefore, higher frequency electroencephalographic activity. Given the hypothesized role of this network in attention 41,42and in “binding” information from multiple cortical areas, 37it is conceivable that loss of synchrony in this network could be associated with sedation, amnesia, and decreased levels of consciousness. In the model of the reticular nucleus of the thalamus, which is known to be associated with sleep and barbiturate anesthesia, 43,47oscillations become more  organized as GABA^Achannels remain open longer. Smaller decrements in the closing rate of the GABA^Achannels are necessary to desynchronize the fast-spiking interneuronal network, whereas larger decrements are necessary to achieve complete synchrony in the model of the reticular nucleus of the thalamus. Thus, different components of the central nervous system may contribute to the overall state, and these contributions, even if involving the same receptor system, may appear at different anesthetic concentrations.

It may be important to appreciate that all aspects of anesthetic action on ion channel behavior and the subsequent impact on individual neurons and the networks they comprise may not be directly related to aspects of the anesthetic state that are clinically recognized and are of direct clinical utility. There is a growing appreciation for the dynamic qualities of the central nervous system, where neurons are continually adjusting their response in accordance with the prevailing pattern of pre- and postsynaptic activity. In various systems, this manifests itself in the form of synaptic plasticity, 105adjustments of dynamic range, 106and alterations in the conductances of populations of ion channels. 107It is conceivable that general anesthesia could affect these processes by altering the pattern of prevailing neural activity and through the dynamic adjustments in response to that activity. These changes could underlie some of the longer-term effects that have been reported following general anesthesia, 108and represent additional features of systems level behavior that may be amenable to computational modeling.

In summary, computational approaches offer additional tools for considering the integrated neural response to anesthetic action, and have demonstrated how anesthetic modulation of ion channel activity could lead to more complex systems level behavior. Collectively, the models described here have exhibited behavior consistent with both excitation and inhibition once anesthetic concentration reaches a clinically relevant threshold. Future progress in this area will require the development of more detailed biophysical models that more realistically incorporate the action of anesthetics and other drugs that modulate anesthetic action. Apart from elucidating the process by which one or more actions at the ion channel lead to general anesthesia, these approaches may eventually pave the way for more rational drug design and use.

The authors thank Marcel Poisit (Research Assistant, Krieger School of Arts and Sciences, Johns Hopkins University, Baltimore, MD) for his assistance with some of the preliminary simulations, and Drs. Roderick G. Eckenhoff (Associate Professor of Anesthesiology, University of Pennsylvania, Philadelphia, PA), Jonas Johansson (Assistant Professor of Anesthesiology, University of Pennsylvania, Philadelphia, PA), and Eric J. Moody (Associate Professor of Anesthesiology and Critical Care Medicine, Johns Hopkins Medical Institutions, Baltimore, MD) for their many suggestions.

Morris-Lecar model of the barnacle giant muscle fiber modified to permit volatile anesthetic modulation of Ca²⁺conductance. 49,56 

TABLE 

C = 20 μF/cm², λ̄^m= 1.0, λ̄^m= 0.1

Initial Conditions :

V(0) =−20 mV, m(0) = 0.065, n(0) = 0.002

Definition of Model Parameters :

I = applied current (μA/cm²)

ḡ^Ca, ḡ^K, ḡ^L= maximal conductance for Ca²⁺, K⁺, and leak currents (mS/cm²)

V = membrane potential (mV)

V^Ca, V^K, V^L= equilibrium potentials for Ca²⁺, K⁺, and leak currents (mV)

m, n = fraction of open Ca²⁺and K⁺channels

m^∞(V), n^∞(V) = fraction of open Ca²⁺and K⁺channels in steady state for given membrane potential

λ^m(V), λⁿ(V) = rate constant for opening of Ca²⁺and K⁺channels (s⁻¹)

λ̄^m, λ̄ⁿ= maximal rate constant for opening of Ca²⁺and K⁺channels (s⁻¹)

V¹, V³= membrane potential at which m^∞(V) and n^∞(V) are equal to 0.5 (mV)

V², V⁴= reciprocal slope of m^∞(V) and n^∞(V) (mV)

C = membrane capacitance (μF/cm²)

a = anesthetic concentration (mM)

a^0.5−Ca= anesthetic concentration at which Ca²⁺-conductance is halved

n^Ca= Hill coefficient for anesthetic modulation of Ca²⁺-conductance

Pinsky-Rinzel Model of hippocampal CA-3 Neuron. 50,51 

TABLE 

Definition of Model Parameters :

V^S, V^D= membrane potentials of somatic and dendritic compartments (mV)

V^Na, V^K, V^Ca, V^L= equilibrium potentials of sodium, potassium, calcium and leak current (mV)

I^S, I^D= applied current to somatic and dendritic compartments (μA/cm²)

I^Na, I^K−DR, I^K−C, I^K−AHP, I^Ca, I^L= sodium, potassium (delayed rectifier), calcium dependent potassium, potassium afterhyperpolarization, calcium, and leak currents (μA/cm²)

m, h = gating variables for activation (m) and inactivation (h) of I^Na

n, s, c, q = gating variables for activation of I^K−DR, I^Ca, I^K−C, I^K−AHP

α^m, β^m, α^h, β^h, αⁿ, βⁿ, α^s, β^s, α^c, β^c, α^q, β^q= kinetic parameters governing activity of the gating variables

m^∞= equilibrium value for gating variable m obtained per Eq. B.9c

χ= saturating function of calcium concentration

ḡ^Na, ḡ^Ca, ḡ^K−DR, ḡ^K−C, ḡ^K−AHP, ḡ^L= maximal conductance of sodium, calcium, delayed rectifier, calcium activated potassium, potassium afterhyperpolarization, and leak currents (mS/cm²)

C = membrane capacitance (μF/cm²)

g^C= coupling strength between somatic and dendritic compartment (mS/cm²)

ρ= fraction of total neuron area in somatic compartment

Golomb-Rinzel model of reticular nucleus of thalamus. 53 

Model Parameters :

V^Ca= 120 mV, V^L=−60 mV, V^syn=−80 mV, ḡ^Ca= 0.5 mS/cm², ḡ^L= 0.05 mS/cm², ḡ^syn= 0.38 mS/cm², α^syn= 1 ms⁻¹, β^syn= .08 ms⁻¹, C = 1 μF/cm², φ= 2 ms⁻¹, N = 100, y = .8

Initial Conditions :

Vⁱ(0) randomly and uniformly distributed between −90 mV and −50 mV (hyperpolarized initial conditions) or −50 mV and −10 mV (depolarized initial conditions). sⁱ(0) and hⁱ(0) obtained from steady state solutions of Eqs. C.2 and C.3for the corresponding value of Vⁱ(0).

Definition of Model Parameters :

Vⁱ= membrane potential of i-th neuron (mV)

N = number of neurons in network

γ= fraction of total number of neurons connected to an individual neuron

S^ij= connection matrix of either zeros or ones, with zeros on the diagonal, where connections were generated randomly so that, on average, γN connections were made to each neuron

V^Ca, V^L, V^syn= equilibrium potential of calcium, leak, and synaptic currents (mV)

ḡ^Ca, ḡ^L, ḡ^syn= maximal conductance of calcium, leak, and synaptic currents (mS/cm²)

hⁱ, sⁱ= gating variable for calcium and synaptic currents for the i-th neuron, representing the fraction of channels that are open

α^syn, β^syn= forward and backward rate constants controlling dynamics of gating variable sⁱwhich represents the fraction of open GABA^Achannels (ms⁻¹)

φ= parameter governing time constant of gating variable for calcium channel

C = membrane capacitance (μF/cm²)

Wang-Buzsáki model of hippocampal interneuronal network. 54 

Model Parameters :

V^Na= 55 mV, V^K=−90 mV, V^L=−65 mV, V^syn=−75 mV, ḡ^Na= 35 mS/cm², ḡ^K= 9 mS/cm², ḡ^L= 0.1 mS/cm², ḡ^syn= 0.1 mS/cm², φ= 5, Θ^syn= 0 mV, C = 1 μF/cm², α^syn= 12 ms⁻¹, β^syn= 0.08 ms⁻¹, γ= 0.6, N = 100, Iⁱ= injection current to i-th neuron (chosen from a Gaussian distribution with mean 1.0 μA/cm²and standard deviation .03 μA/cm²)

Initial Conditions :

Vⁱ(0) randomly and uniformly distributed between −70 mV and −50 mV. hⁱ(0), nⁱ(0), and sⁱ(0) obtained from steady state solutions of Eqs. D.2–4for the corresponding value of Vⁱ(0).

Definition of Model Parameters :

Vⁱ= membrane potential of i-th neuron (mV)

N = number of neurons in network

γ= fraction of total number of neurons connected to an individual neuron

S^ij= connection matrix of either zeros or ones, with zeros on the diagonal, where connections were generated randomly so that, on average, γN connections were made to each neuron

V^Na, V^K, V^L, V^syn= equilibrium potential of sodium, potassium, leak, and synaptic currents (mV)

ḡ^Na, ḡ^K, ḡ^L, ḡ^syn= maximal conductance of sodium, potassium, leak, and synaptic currents (mS/cm²)

mⁱ, hⁱ, nⁱ, sⁱ= gating variables governing activation and inactivation of sodium current, activation of potassium current, and activation of inhibitory GABA^Acurrent

m^∞= steady state value of gating variable governing activation of sodium current

φ= dimensionless parameter governing time constant of gating variable governing sodium inactivation and potassium activation

α^h, β^h, αⁿ, βⁿ= kinetic parameters governing activity of the gating variables hⁱ, and nⁱ

α^syn, β^syn= forward and backward rate constants controlling dynamics of gating variable sⁱwhich represents the fraction of open GABA^Achannels (ms⁻¹)

Iⁱ= injection current to i-th neuron (μA/cm²). Injection currents, although constant throughout the simulation for each neuron, were not uniform and were chosen from a Gaussian distribution.

F(V^j) = modulates synaptic opening as a function of presynaptic membrane potential

Θ= presynaptic membrane potential at which synaptic opening rate reaches its half-maximal value

C = membrane capacitance (μF/cm²)