A Mechanistic Neural Field Theory of How Anesthesia Suppresses Consciousness: Synaptic Drive Dynamics, Bifurcations, Attractors, and Partial State Equipartitioning

With the advances in biochemistry, molecular biology, and neurochemistry there has been impressive progress in understanding the molecular properties of anesthetic agents. However, there has been little focus on how the molecular properties of anesthetic agents lead to the observed macroscopic property that defines the anesthetic state, that is, lack of responsiveness to noxious stimuli. In this paper, we use dynamical system theory to develop a mechanistic mean field model for neural activity to study the abrupt transition from consciousness to unconsciousness as the concentration of the anesthetic agent increases. The proposed synaptic drive firing-rate model predicts the conscious–unconscious transition as the applied anesthetic concentration increases, where excitatory neural activity is characterized by a Poincaré–Andronov–Hopf bifurcation with the awake state transitioning to a stable limit cycle and then subsequently to an asymptotically stable unconscious equilibrium state. Furthermore, we address the more general question of synchronization and partial state equipartitioning of neural activity without mean field assumptions. This is done by focusing on a postulated subset of inhibitory neurons that are not themselves connected to other inhibitory neurons. Finally, several numerical experiments are presented to illustrate the different aspects of the proposed theory.


Introduction
The in vitro effects of anesthetic agents have been intensely investigated, and the actions of anesthetic agents on single neurons have been described [1,2]. There is compelling evidence that at least some modern anesthetics alter postsynaptic potentials. A key goal of anesthetic and neuroscience research is to understand how this effect on the single neuron translates into an abrupt transition from consciousness to unconsciousness as the concentration of the agent increases [3][4][5][6][7].
A rigorous analysis of the effect of anesthetic agents ideally would be within the context of the Hodgkins-Huxley nonlinear differential equations describing the electrical characteristics of neurons [8,9]. However, the complexity of this model limits its tractability for systems of multiple neurons. Since the seminal publication of Wilson and Cowan [10] it has been common to instead utilize firing-rate models. However, even with the considerable mathematical simplifications afforded by the assumptions of firing-rate models, the immense dimensionality of the brain requires further assumptions and approximations. Typically, firing-rate models are simplified further by either mean field assumptions, which postulate that the brain is organized into a limited number of populations of identical neurons, or by assuming that the strength of connections between neurons have a specific (typically normal) distribution. In some instances the firing-rate model is further simplified by postulating specific network architectures [11].
While firing-rate models have been extensively investigated in the neuroscience literature, there has been limited applications of these models to the effects of anesthetic agents. One of the earlier investigations came from Steyn-Ross et al. [12], who extended a generalized mean field firing-rate model that postulates specific short range and long range connections between cortical neurons. This model was originally developed by Liley et al. [13] to analyze the effects of anesthetic agents on the mean electrical potential of the cell body of neurons. The model proposed in [12] predicts that the transition from consciousness to unconsciousness induced by anesthetic agents can be characterized as a first-order phase transition to a state of decreased neural activity. However, the significance of decreased neural activity for the transition from consciousness to unconsciousness is unclear.
Anesthetic agents typically cause paradoxical increases in the activity of the brain, sometimes to the point of seizure activity, before leading to decreased activity. A primary experimental modality for investigating the induction of anesthesia is the electroencephalogram (EEG). This is manifested in the EEG where anesthetic agents typically cause an increase in the amplitudes of certain frequency bands at the point of transition to unconsciousness with subsequent burst suppression and eventually a flat line EEG as the anesthetic concentration increases. There has been significant interest in understanding the effects of anesthetic agents on the EEG [2][3][4][14][15][16][17] and the theoretical models can predict this paradoxical increase in neural activity observed experimentally.
In contrast to the phase transition model proposed in [12], another general theory of the loss of consciousness is that it reflects loss of information flow between different parts of the brain [18]. John and Prichep [19] have proposed a model of specific sequential neuroanatomic targets for anesthetic agents, originating in the reticular activating system and progressing through the mesolimbic system, with closure of thalamic gates and blockade of thalamocortical reverberations and uncoupling of parietal-frontal interactions with loss of consciousness. It is plausible that a decrease of neural activity in specific segments or neuronal populations of the brain could lead to the loss of information flow and there is experimental evidence to support this hypothesis [20].
In this paper, we use dynamical system theory to develop a mechanistic mean field model for neural activity. Our focus here is to understand general neural activity without attempting to predict the anatomic specificity of the John and Prichep theory [19]. Specifically, using the notion of synaptic drive-a measure of neuronal population activity that captures the influence of a given neuron population on the behavior of the neuronal network-and topological system stability, we show that this simplified model predicts decreasing mean excitatory neural activity with increasing anesthetic potency at the level of the single neuron. We seek conditions that lead to decreasing excitatory neural activity despite the reservations noted above around the role of decreasing activity for the induction of anesthesia for two reasons. First, as noted above, it is plausible that decreasing excitatory activity in selected subpopulations of the brain could lead to the interruption of information flow that may characterize unconsciousness. Second, anesthetic agents do eventually lead to decreased activity as the concentration of the anesthetic agent increases. While we would expect this of any phenomenological model for addressing how anesthesia suppresses consciousness, the novelty of the model is the additional prediction of a conscious-unconscious transition as the applied anesthetic concentration increases, where excitatory neural activity is characterized by a stable limit cycle. Then, as the anesthetic concentration increases further, the proposed dynamical system model undergoes a supercritical reverse Hopf bifurcation and transitions from a limit cycle behavior to an asymptotically stable equilibrium corresponding to an anesthetized state, where the excitatory neuron activity is very low.
Next, we extend our mean field theory to incorporate a more complex model for the postsynaptic potential. Specifically, by using an Erlang-type time multiplied exponential decay model for the postsynaptic potential rather than a simple exponential decay [11,13], we can account for the delay in peak amplitude of the postsynaptic potential that occurs after a neuron discharges ( [11], Fig. 8.5). In this case, we see biphasic responses in the mean neural activity that seems to parallel the rise and decay of the postsynaptic potential. The mean field synaptic drive model we present is based on the assumption that, within the populations of excitatory and inhibitory neurons, second-order terms reflecting the variation of synaptic connection strengths and the variation of synaptic drives can be ignored [5].
Finally, we investigate the synaptic drive model without mean field assumptions by postulating the existence of a subpopulation of inhibitory neurons that themselves do not receive inputs from other inhibitory neurons. In this case, we provide sufficient conditions for global asymptotic and exponential partial synaptic drive equipartitioning for our excitatory and inhibitory cortical neuronal network. Specifically, we show that as the inhibitory time constants increase (one of the demonstrated in vitro effects of some anesthetic agents [1,2]), the excitatory neurons that are coupled to inhibitory neurons all approach a zero synaptic drive state, whereas the inhibitory neurons that themselves are not coupled to inhibitory neurons converge to a constant synaptic drive state.
The notation used in this paper is fairly standard. Specifically, for x ∈ R n , we write x ≥≥ 0 (respectively, x 0) to indicate that every component x i of x is nonnegative (respectively, positive). In this case, we say that x is nonnegative or positive, respectively. Likewise, A ∈ R n×m is nonnegative or positive if every entry of A is nonnegative or positive, respectively, which is written as A ≥≥ 0 or A 0, respectively. For A ∈ R n×n , we write A ≥ 0 (respectively, A > 0) to indicate that A is nonnegative (respectively, positive) definite. Furthermore, we write R n + and R n + to denote the nonnegative and positive orthants of R n , that is, if x ∈ R n , then x ∈ R n + and x ∈ R n + are equivalent, respectively, to x ≥≥ 0 and x 0. In addition, we write 0 n×m to denote the n × m zero matrix, I n to denote the n × n identity matrix, and e n ∈ R n to denote the ones vector of order n, that is, e n = [1, 1, . . . , 1] T ; if the order of e n is clear from the context we simply write e for e n . Finally, we write (·) T to denote the transpose, tr(·) to denote the trace, det(·) to denote the determinant, and (·) to denote the Fréchet derivative.

A Model for Excitatory and Inhibitory Neural Populations
Biological neural network models predict a voltage in the receiving or postsynaptic neuron given by [7,9] v X where v X i (·), i = 1, . . . , n X , X ∈ {E, I}, is the excitatory (X = E) and inhibitory (X = I) voltage in the ith receiving neuron, A XY ij , X, Y ∈ {E, I}, are constants representing the coupling strengths (in volts) of the j th neuron on the ith neuron, k, k = 1, . . . , enumerate the action potential or firings of the excitatory and inhibitory transmitting (presynaptic) neurons at firing times t k and t k , respectively, and α E j (·) and α I j (·) are dimensionless functions describing the evolution of the excitatory and inhibitory postsynaptic potentials, respectively. Using a (possibly discontinuous) function f i (·) to represent the firing rate (in Hz) of the ith neuron that determines the instantaneous number of spikes that a neuron fires over an infinitesimal time interval [t, t + dt], and assuming that the firing rate is a function of the voltage v E i (·) (resp., v I i (·)) across the membrane of the ith neuron given by where f j (v X j (t)) dt, X ∈ {E, I}, is the probability of a spike occurring in the time interval [t, t + dt], and A XY ij , X, Y ∈ {E, I}, are neural connectivity weights, with units of volts, representing the coupling strength of the j th neuron on the ith neuron such that A XE ij > 0 and A XI ij < 0, X ∈ {E, I}, if the j th neuron is connected (i.e., contributes a postsynaptic potential) to the ith neuron, and A XY ij = 0, otherwise. Alternatively, f i (·) can represent an ensemble average of firing rates of neurons in a single population i over [t, t + dt]. In this case, the firing rate of the population is exactly the same as that of an individual neuron and (2) and (3) represent excitatory and inhibitory voltages for the ith population. Furthermore, v E thi (·) and v I thi (·) are continuous input threshold voltages characterizing nerve impulses from sensory (pain) receptors, thermo (temperature sensing) receptors, or proprioceptive (motion sensing) receptors. Alternatively, v E thi (·) and v I thi (·) can be thought of as inputs from the reticular activating system within the brainstem responsible for regulating arousal and sleep-wake transitions. Note that A EE ii A II ii 0 by definition. Next, defining the synaptic drive-a dimensionless quantity-of each (excitatory or inhibitory) neuron as the convolution of the presynaptic firing rate with the postsynaptic potential given by [7,9] and assuming an exponential decay of the synaptic voltages of the form [9] α (E,I) i  (4) and (5) that A more general Erlang-type time multiplied exponential decay model for the postsynaptic potential is considered in Sect. 5. Note that the synaptic drive accounts for pre-and postsynaptic potentials in a neural network to give a measure of neural activity in the network. In particular, it follows from (4) that the synaptic drive quantifies the present activity (via the firing rate) along with all previous activity within a neural population appropriately scaled (via a temporal decay) by when a particular firing event occurred. Hence, the synaptic drive provides a measure of neuronal population activity that captures the influence of a given neuron population on the behavior of the network from the infinite past to the current time instant.
The above analysis reveals that a form for capturing the neuroelectric behavior of biological excitatory and inhibitory neuronal networks can be written as where S i (t) ∈ D ⊆ R, t ≥ 0, is the ith synaptic drive, v thi (t) ∈ R, t ≥ 0, denotes the input threshold voltage to the ith neuron, A ij is a constant representing the coupling strength of the j th neuron on the ith neuron, τ i = λ i is a synaptic time scale, B i is a constant gain for the firing rate of the ith neuron, and f i (·) is a nonlinear activation function describing the relationship between the synaptic drive and the firing rate of the ith neuron. Even though the computation of the neural connectivity matrix and the graph topology of human brain networks are intractable, the phenomenological firing-rate model (6) and (7) involving the averaged behavior of the spiking rates of groups of neurons can be used to predict network system changes due to changes in a subset of the neuronal connectivity matrix A XY containing the entries A XY ij in order to understand how the large neuron population changes qualitatively with the induction of anesthesia. The relevance of the model to realistic systems can be appraised by its prediction of salient aspects of anesthesia. In particular, predicting the abrupt transition from consciousness to unconsciousness resembling a phase transition as well as the biphasic response during induction of anesthesia leading to a paradoxical phase prior to the loss of consciousness [21].
In such population models, the activity of a neuron (population), that is, the rate at which the neuron (population) generates an action potential (i.e., "fires") is modeled as a function of the voltage (across the membrane). In this paper, we will assume continuous half-wave rectification activation functions as well as smooth sigmoidal functions. Specifically, for a typical neuron [22] where i ∈ {1, . . . , n} and [x] + = x if x ≥ 0, and [x] + = 0, otherwise. The activation function (9) reflects the fact that as the voltage increases across the membrane of the ith neuron, the firing rate increases as well. Often, the membrane potential firingrate curve exhibits a linear characteristic for a given range of the voltages. At higher voltages, however, a saturation phenomenon appears, indicating that the full effect of the firing rate has been reached. To capture this effect, f i (·) can be modeled as a smooth (i.e., infinitely differentiable) sigmoidal function where i ∈ {1, . . . , n}, γ 0, and f max = lim x→∞ f i (x) denotes the maximum firing rate.

A Two-Class Mean Excitatory and Mean Inhibitory Synaptic Drive Model
To avoid the complexity of the large-scale neural network model (6) and (7), in this section we consider a mean field model. Specifically, the excitatory and inhibitory synaptic drive model given by (6) and (7) can be reduced to a two-class mean excitatory and mean inhibitory model. In particular, with continuously differentiable (6) and (7) collapse to (see [7] for details) where a j =1 S I j (t) denote the mean excitatory synaptic drive and mean inhibitory synaptic drive in dimensionless units, respectively. Note that the constants a, b, c, and d are nonnegative. Equations (11) and (12) represent the spatial average (mean) dynamics of the system given by (6) and (7), and they are predicated on a mean field assumption that reduces the complex (approximately 10 11 × 10 11 ) neuronal connectivity matrix to a 2 × 2 excitatory-inhibitory system. This is a drastic assumption, but one which has been commonly used in theoretical neuroscience going back to the pioneering work of Wilson and Cowan [10].
To study the dynamic behavior of the mean field model (11) and (12), we assume the postsynaptic activation function is given by (10). The next set of propositions and theorem present several results on the dynamic behavior of (11) and (12). Here, we use the language of topological dynamics (i.e., flows, equilibria, periodic orbits, limit sets) to study the dynamic behavior of (11) and (12). Recall that a set D ⊆ R 2 is positively invariant with respect to (11) and (12) if D contains the orbits of all its points. Moreover, recall the standard Lyapunov and asymptotic stability definitions for an equilibrium point of (11) and (12) given in [23], and recall that if all solutions of (11) and (12) are bounded, then it follows from the Peano-Cauchy theorem ( [23], p. 76), that the maximal solution to (11) and (12) exists on the semi-infinite interval [0, ∞), and hence, (11) and (12) are forward complete.
The following propositions are needed for the main result of this section.
Proposition 1 [7] Consider the two-class mean excitatory and mean inhibitory synaptic drive network given by (11) and (12). If S Hence, the vector field f along the line S I = 0 is directed inward toward the region M. Finally, let [0 −1] T be a normal vector along the line S I = f max λ I directed toward the region M and note that, Hence, the vector field f along the line S I = f max λ I is directed inward toward the region M. Thus, along the boundary ∂M of M the vector field f is directed inward toward the region M, and hence, M is positively invariant with respect to (11) and (12).
A visualization of the region M is shown in Fig. 1. Note that the equilibrium points (S E e , S I e ) of (11) and (12) are characterized by the solution to f cS Since f (·) is increasing, f (·) is invertible. Hence, it follows from (15) and (16) that  Proposition 3 Consider the two-class mean excitatory and mean inhibitory synaptic drive network given by (11) and (12).
then there exist input voltages v I th and v E th such that (11) and ( Note that the maximum values of the gradients of g(·) and h −1 (·) occur at the inflection points of g(·) and h(·), respectively, and from (19) and (20) the inflection points of g(·) and h(·) correspond to S E = f max λ E 2 and S I = f max λ I 2 , respectively. Thus, the maximum values of the gradients of g(·) and h −1 (·) correspond to 1 b (a − 4 γf max λ E ) and , respectively. Now, the condition that the maximum values of the gradient of h −1 (·) is greater than the maximum value of the gradient of g(·) implies that Hence, there exist input voltages v I th and v E th such that (11) and (12)  Theorem 1 Consider the two-class mean excitatory and mean inhibitory synaptic drive network given by (11) and (12), and assume that the hypothesis in Proposition 3 is satisfied. If tr J (S E e , S I e ) < 0, then there exist input voltages v I th and v E th such that (11) and (12) (11) and (12) possess a limit cycle. Moreover, the limit cycle is stable.
Proof It follows from Proposition 3 that (11) and (12)  is an unstable equilibrium point, it follows that R is positively invariant with respect to (11) and (12). Hence, it follows from the Poincaré-Bendixson theorem ( [23], p. 109) that R contains a periodic orbit. Furthermore, since R is positively invariant and the equilibrium point (S E e , S I e ) is unstable, the limit cycle is stable.

Large Neural Populations, Synchronization, and Partial State Equipartitioning
One of the most important questions in neuroscience is how do neurons, or collections of neurons, communicate. There is extensive experimental verification that collections of neurons may function as oscillators and the synchronization of oscillators may play a key role in the transmission of information within the central nervous system [24][25][26]. This may be particularly relevant to understanding the mechanism of action for general anesthesia. It has been known for a long time that general anesthesia has profound effects on the spectrum of oscillations in the electroencephalograph (EEG) [15,27,28]. More recently, the authors in [19] have suggested that thalamocortical circuits function as neural pacemakers and that alterations in the thalamic oscillations are associated with the induction of general anesthesia. Furthermore, it is well known that anesthetic drugs frequently induce epileptiform activity as part of the progression to the state of unconsciousness [29]. Multiple lines of evidence indicate that anesthetic agents impact neural oscillators. In addition, epileptiform activity implies synchronization of oscillators. This leads to the possibility that synchronization of these oscillators is involved in the transition to the anesthetic state.
Sufficient conditions for synchronization and full state equipartition for an excitatory and inhibitory cortical neural network are given in [7]. It has been observed, however, that when patients lose consciousness there are differences in neural activity in different anatomic regions of the cerebral cortex. This can be captured by biological neural network models that exhibit partial synchronization, wherein part of the system's state (neural activity) is synchronized and the other parts fire at normal levels.
Before proceeding, it is important to distinguish between the notions of synchronization and state equipartitioning or consensus. System synchronization refers to the fact that the dynamical system states achieve temporal coincidence over a finite or infinite time horizon, whereas state equipartitioning refers to the fact that the dynamical system states converge to a common value over a finite or infinite time horizon. Hence, both notions involve state agreement in some sense. However, equipartitioning involves convergence of the state values to a constant state, whereas synchronization involves agreement over time instants. Thus, equipartitioning implies synchronization; however, the converse is not necessarily true. It is only true in so far as consensus is interpreted to hold over time instants.
In order to address partial synchronization, we make a high-level assumption regarding neuronal connectivity by postulating the existence of a subset of inhibitory neurons that themselves do not receive inhibitory input. It is important to note that the existence of such a subset has not been demonstrated experimentally. However, it is not biologically implausible and, as we shall demonstrate, this assumption leads to sufficient conditions for partial synchronization.
To address the problem of partial synchronization, let S E (t) . . , n I , and assume B E i = B I i = 1 so that the vector-matrix form of (6) and (7) can be written aṡ where denote the vector activation functions describing the relationship between the synaptic drives and the firing rates of neurons with f i (x), x ∈ R, i = 1, . . . , n E , and f k (x), x ∈ R, k = 1, . . . , n I , defined as in (9).
. . , S I q (t)] T denotes the vector of inhibitory synaptic drives that are coupled to both excitatory and inhibitory neurons andS I 2 (t) [S I q+1 (t), . . . , S I n I (t)] T denotes the vector of inhibitory synaptic drives that are only coupled to excitatory neurons. In this case, the neuronal connectivity matrix A can be partitioned as where (25) and (26) can be written aṡ where (31), and (32) can be written asṠ whereÃ The following proposition and definitions are needed for the statement of the main results of this section.

Definition 2
The biological neural network given by (34) and (35) is said to be globally asymptotically partially synchronized if

Definition 3
The biological neural network given by (34) and (35) is said to be globally exponentially partially synchronized if there exist constants p > 0 and ρ > 0 such that Remark 1 It follows from Definitions 1 and 3 that if the biological neural network given by (34) and (35) is globally exponentially stable with respect to S uniformly iñ S I 20 , then (34) and (35) is globally exponentially partially synchronized.
For the statement of the next theorem, N (X) denotes the nullspace of the matrix X. (34) and (35) with the vector activation functions defined by (27) and (28).
Proof Consider the partial Lyapunov function candidate V : R n E +q + → R given by where λ max (P ) and λ min (P ) denote the maximum and minimum eigenvalue of P , respectively, and · 2 denotes Euclidean norm. It follows that the derivative of V (S) along the trajectories of (34) is given bẏ Next, it follows from (36) that or, equivalently, Now, using (40) it follows from (38) thaṫ or, equivalently, Now, consider the dynamical systeṁ and note that Using the comparison principle ( [23], p. 126), it follows that . . , n I , or, equivalently, Next, it follows from (41) thaṫ Since . . , n E + q, given by (9), and R ∈ R (n E +q)×(n E +q) is a diagonal positive-definite matrix, it follows that Now, it follows from (37) and (47) thaṫ Alternatively, if Ω ≥ 0 and N (Ω) = span(e n E +q ) holds, thenV (S(t),S I 2 (t)) ≤ 0, t ≥ 0, and hence, V (S(t)) ≤ V (S 0 ) for all t ≥ 0. Next, since P is positive definite anḋ V (S(t),S I 2 (t)) is a nonincreasing function of time, it follows that V (S(t)) is bounded for all t ≥ 0, and hence, S(t) is bounded for all t ≥ 0, which further implies that is bounded for all t ≥ 0, and hence, there existsf I 2 max such thatf Next, consider the dynamical systeṁ and note that since L I 2 is a diagonal positive-definite matrix, it follows that x(t) is bounded for all t ≥ 0. Using the vector comparison principle [23], it follows from (50) and (51) thatS I 2 (t) ≤≤ x(t), t ≥ 0, and hence,S I 2 (t) is bounded for all t ≥ 0. Since S(t) andS I 2 (t) are bounded for all t ≥ 0, it follows thatV (S(t),S I 2 (t)) is bounded for all t ≥ 0. Thus,V (S(t),S I 2 (t)) is uniformly continuous in t. Finally, it follows from Barbalat's lemma ( [23], p. 221) thatV (S(t),S I 2 (t)) → 0 as t → ∞, which, since N (Ω) = span(e n E +q ), implies that (34) and (35) is globally asymptotically partially synchronized.
There is a body of experimental evidence that the mechanism of action of some anesthetic drugs is the prolongation of the inhibitory postsynaptic potential time constant [1,2]. While Theorem 2 provides sufficient conditions for partial synchronization, it is not readily apparent how these conditions may be met as the inhibitory postsynaptic potential time constant increases. Next, we address this question for the specific activation function [7] where f max denotes the maximum firing rate.
Condition (53) is satisfied if A EI 2 does not contain any row of zeros and L I 2 is small enough or, equivalently, the time constants of the inhibitory neurons that themselves do not receive any inhibitory inputs are large enough. This condition captures an important aspect of the effect of the anesthetic cascade. In particular, as the anesthetic drug concentration increases, resulting in an increase in the inhibitory time constants, the excitatory synaptic drives that receive inhibitory inputs from inhibitory neurons, which themselves do not receive inhibitory inputs, converge to zero, reflecting a state of unconsciousness. Alternatively, if only a subset E 1 of excitatory neurons receive inhibitory inputs from the inhibitory neurons in I 1 , then it can be shown that only the synaptic drives of the excitatory neurons in E 1 will converge to zero as the time constants of the inhibitory neurons in I 1 increase.

Generalized Neural Population Models, Nonmonotonic Postsynaptic Potentials, and Drug Biphasic Response
The administration of increasing anesthetic doses can lead to a paradoxical state of excitement in the patient prior to decreases in the level of consciousness. This paradoxical boost in brain activity prior to hypnotic induction is known as drug biphasic response [21]. There is also a second biphasic surge in the EEG power as the patient emerges from unconsciousness. The model proposed in Sect. 2 does not capture such phenomena. Models that predict the aforementioned characteristics are of great clinical importance in providing the phenomenological trends of the anesthetic cascade.
In this section, we incorporate a more complex model for the postsynaptic potential to account for the delay in peak amplitude of the postsynaptic potential that occurs after a neuron discharges ( [11], Fig. 8.5).
To capture an initial overshoot in the synaptic drive dynamics in an attempt to see biphasic responses in our neural field model, we assume the functional form of the postsynaptic potential is given by In this case, it follows from (4) and (63) that Now, defining it follows that Next, it follows from (64) and (65) that Differentiating (67) with respect to time, we obtain Now, using the expressions for the excitatory and inhibitory voltage given by (2) and (3), respectively, it follows that The above analysis reveals that an alternative form to (8) for capturing the neuroelectric behavior of biological excitatory and inhibitory neuronal networks can be given by the second-order differential equations To address the problem of partial synchronization for the model given by (69) and (70), the vector-matrix form of (69) and (70) can be written as where all vector-matrix notation are defined as in Sect. 4. Once again, we partition the state vector S I (t), t ≥ 0, as in Sect. 4 so that the neuronal connectivity matrix is partitioned as in (29), and hence, (72) and (73) can be written as or, equivalently,
Condition (98) implies that all the excitatory neurons receive inhibitory inputs from the inhibitory neurons that themselves do not receive inhibitory inputs. As in the case of the simpler neural population model discussed in Sect. 4, if only a subset of excitatory neurons receive the inhibitory inputs from the inhibitory neurons that themselves do not receive inhibitory inputs, then only the synaptic drives of this subset of excitatory neurons will converge to zero as the inhibitory time constants increase.

A Second-Order Mean Excitatory and Inhibitory Synaptic Drive Model
In this section, we reduce the model given in (69) and (70) to a mean excitatory and mean inhibitory model. Consider (69) and (70) with continuously differentiable (115) Using the average and perturbed expression for A XY ij , X, Y ∈ {E, I}, (113) and (114) can be written as where S denote the mean excitatory synaptic drive and mean inhibitory synaptic drive in dimensionless units, respectively.
where δ E i (t) and δ I i (t) are deviations from the mean, (116) and (117) become Next, assume that all terms with a factor XY ij , X, Y ∈ {E, I}, i = 1, . . . , n X and j = 1, . . . , n Y , in (118) and (119) are small relative to the remaining terms in f (·). Then a first-order expansion of (118) and (119) gives Now, assuming that the higher-order terms can be ignored, (120) and (121) become Finally, summing (122) and (123) over i = 1, . . . , n E and i = 1, . . . , n I , dividing by n E , and n I , respectively, using (115), and assuming v E thi (t) = v E th and v I thi (t) = v I th , t ≥ 0, it follows that the average excitatory synaptic drive and the average inhibitory synaptic drive are given by A similar analysis to the analysis provided in Sect. 3 can be carried out for the coupled second-order differential equations (124) and (125). This can allow us to qualitatively study the abrupt transition from consciousness to unconsciousness as well as the biphasic response during induction of anesthesia. Initial quantitative numerical experiments of this model are presented in Sect. 7.

Illustrative Numerical Simulations of the Neural Field Model
In this section, we use the different analysis results developed in the paper to present several numerical experiments that illustrate the qualitative behavior of the proposed firing-rate models.

Two-Class Mean Excitatory and Mean Inhibitory Synaptic Drive Model
First, we consider the two-class mean field model (11) and (12) with a single equilibrium point. Figures 3 and 4 show the case when the equilibrium point is unstable, and hence, the system possesses a stable limit cycle. Next, we illustrate the effect of the excitatory sensory input v E th on the mean excitatory and mean inhibitory synaptic drives. Here, we fix the parameters to be as in the simulation shown in Fig. 3 and vary v E th . It can be seen from Fig. 5 that, for the value of v E th below v E th1 = −1.6, (11) and (12) have exactly one asymptotically stable equilibrium point. As v E th increases to v E th2 = −1.6, the equilibrium point becomes unstable. At this point, a supercritical Hopf bifurcation occurs, giving rise to a stable limit cycle. As v E th increases beyond v E th2 , the unstable equilibrium point reverts to an asymptotically stable equilibrium point and the stable limit cycle disappears. The value of S E e and S I e increases as v E th increases; see Fig. 5. The effect of v E th on the mean firing rates and mean synaptic drives of the excitatory and the inhibitory neurons is shown in Figs. 6 and 7, respectively. As v E th increases, the mean firing rates and mean synaptic drives for the excitatory and the inhibitory neurons increase for the chosen parameter values.  (11) and (12) for a = 10, b = 9, The input voltage v I th has a negative effect on the mean excitatory firing rate and synaptic drive. As can be seen from Fig. 8, when v I th increases the mean excitatory firing rate and synaptic drive decrease for the chosen parameter values.
Finally, we relate the effect of an anesthetic drug, corresponding to an increase in the inhibitory time constant λ I , to the mean excitatory and mean inhibitory synaptic drives. Here, we fix the parameters to be as in the simulation shown in Fig. 3 and vary λ I , which corresponds to concentration variations in the anesthetic drug. As can be seen in Fig. 9, for the value of λ I below λ I 1 = 0.85, (11) and (12) have exactly one asymptotically stable equilibrium point. As λ I increases to λ I 1 , the equilibrium point becomes unstable. At this point, a supercritical Hopf bifurcation occurs, giving rise to a stable limit cycle. As λ I increases beyond λ I 2 = 2.3, the unstable equilibrium point reverts to an asymptotically stable equilibrium point and the stable limit cycle disappears. The value of S E e decreases as λ I increases and approaches zero as λ I becomes large; see Fig. 9. where the single asymptotically stable equilibrium point becomes unstable and gives rise to a stable limit cycle. At v E th = 0.6, the unstable equilibrium point reverts to an asymptotically stable equilibrium point and the stable limit cycle disappears The mean excitatory synaptic drive S E (t), t ≥ 0, also decreases as λ I increases and approaches zero with increasing λ I ; see Fig. 10. The effect of anesthetic drug can be explained by observing that if λ I is below λ I 1 , corresponding to no drug effect, then the mean excitatory synaptic drive is at a high level, which corresponds to an awake state. The range of drug concentrations where we cannot predict with certainty whether the patient will respond to noxious stimuli or not could certainly be explained by stochastic variation in v E th and v I th ; the postulated inputs from pain, sensorimotor, and proprioceptive sensors, or the reticular activating system.
However, an alternative explanation is that during the induction of anesthesia, the patient may or may not respond to a certain range of anesthetic drug concentration. This is captured in our model by the existence of limit cycle for a range of λ I between λ I 1 and λ I 2 . In this case, the mean excitatory synaptic drive has a maximum and a minimum value as can be seen in Fig. 10. The probability of observing the response of the patient depends on the time in the cycle at which the observation is made. The patient will respond if the observation occurs at the time when the mean excitatory synaptic drive is at its peak and will not respond if the mean excitatory synaptic drive is at its nadir. If the concentration of the anesthetic drug is further increased, corresponding to an increase in λ I , then the mean excitatory synaptic drive approaches zero as shown in Fig. 10. At this point, the patient is deeply sedated and will not respond.

Partial Synchronization of Synaptic Drive Dynamics
For our first example, we demonstrate partial synchronization for the model given by (34) and (35) with three excitatory neurons E 1 -E 3 and three inhibitory neurons I 1 -I 3  Fig. 11. The neural connectivity matrix A is given by which implies that the two inhibitory neurons I 2 and I 3 do not receive inhibitory inputs. Here, we assume that all the excitatory neurons have the same  We use the MATLAB ® toolbox YALMIP for solving the linear matrix inequalities (LMIs) given by (36)   Hence, the biological neural network given in (34) and (35) is globally exponentially partially synchronized. As can be seen in Fig. 12, the synaptic drive of the excitatory neurons E 1 to E 3 and the inhibitory neuron I 1 converges to zero, whereas the synaptic drive of the inhibitory neurons I 2 and I 3 that themselves do not receive inhibitory inputs do not converge to zero. Next, we increase the number of neurons in the network to twelve with six excitatory neurons E 1 -E 6 and six inhibitory neurons I 1 -I 6 . The connections between neurons are also increased as shown in the neural connectivity matrix A given by Fig. 9 Bifurcation diagram of (11) and (12) for a = 10, b = 9, c = 6, d = 1, v E th = −0.5, v I th = −2.5, λ E = 1, f max = 1, and γ = 1 with λ I as a bifurcation parameter. The solid line represents asymptotically stable equilibrium point, whereas the dashed line represents unstable equilibrium point. A supercritical Hopf bifurcation occurs at λ I 1 = 0.85 where the single asymptotically stable equilibrium point becomes unstable and gives rise to a stable limit cycle. At λ I 2 = 2.3, the unstable equilibrium point reverts to an asymptotically stable equilibrium point and the stable limit cycle disappears which implies that the three inhibitory neurons I 4 , I 5 , and I 6 do not receive inhibitory inputs. Here, we assume that all the excitatory neurons have the same time constant λ E = 0.05 s and all the inhibitory neurons have the same prolonged time constant  Hence, the biological neural network given in (34) and (35) is globally exponentially partially synchronized. As can be seen in Fig. 13, the synaptic drive of the excitatory neurons E 1 to E 6 and three of the inhibitory neurons I 1 , I 2 , and I 3 converges to zero,     Hence, the biological neural network given in (79) and (80) is globally exponentially partially synchronized. As can be seen in Fig. 16, the synaptic drive of the excitatory neurons E 1 to E 4 and two of the inhibitory neurons I 1 and I 2 converges to zero, whereas the synaptic drive of the inhibitory neurons I 3 and I 4 do not converge to zero.
Next, we apply Theorem 5 to the connectivity matrix given by (127 It can also be seen from Figs. 16 and 18 that the synaptic drive of the excitatory neurons E 1 to E 4 increase initially before converging to zero. This potentially captures the drug biphasic response, where there is a boost in brain activity prior to hypnotic induction. Finally, we further explore the drug biphasic phenomena using the second-order mean excitatory and mean inhibitory synaptic drive model given by (124) and (125) with the activation function f (·) given by (10). Once again, we observe that there is an initial boost in the mean excitatory synaptic drive just before convergence to zero; see Fig. 19.

Conclusion and Discussion
Over the past decade there have been significant advances in our understanding of how anesthetic agents affect the properties of neurons [1,2]. This has led, in turn, to efforts to understand how molecular mechanisms translate into the induction of general anesthesia at the macroscopic level [3][4][5][6][7][12][13][14][15][16][17]. In particular, since the pioneering work of Steyn-Ross et al. [12] there has been particular focus on how the molecular properties of anesthetic agents lead to the observed changes in the electroencephalogram that are associated with the induction of anesthesia [12][13][14][15][16][17]. These  1.2, 1.3, 1.25], and f max = 0.5. The synaptic drive of the excitatory neurons E 1 to E 4 do not converge to zero models, of necessity, have focused on the electrochemical potential of the neuron and the analyses consider both the temporal and the spatial domain, specifically invoking a continuous spatial domain to account for interneuron connectivity. The complexity of the brain requires specific assumptions as regards this connectivity.
In this paper, we used a synaptic drive formulation of neuronal activity to derive a mechanistic model that can be employed to understand the effects of anesthetic agents on neuronal firing rates. Our focus is on demonstrating a decrease in the firing rate of excitatory neurons as a consequence of the known effects of anesthetic agents on single neurons. Thus, there is the implicit assumption that anesthesia is a result of decreased neuronal activity. We do note that this may be overly simplistic since the induction of anesthesia is associated with an initial increase in neuronal activity. However, we seek to understand how decreased activity might develop because (i) anesthetic agents at sufficiently high concentrations will cause a decrease in neuronal activity and (ii) one theory of anesthesia is that it represents the interruption of flow of "information" through the brain and this could certainly result from decreased neuronal activity in specific regions.
Even this seemingly simple effect cannot be demonstrated using the full synaptic drive model. Our model begins with a spatially discrete description of neuronal connectivity. However, given the immense dimensionality of the brain, in the first part of the paper we simplify our model using a mean field assumption. Specifically, we assume that within the populations of excitatory and inhibitory neurons, second-order terms that are the product of variances of synaptic connection strengths and variances of synaptic drives (from the mean) can be ignored [5].
We show that even this simplified two-state mean field theory exhibits interesting behavior. In particular, for certain values of model parameters, the mean excitatory synaptic drive has an asymptotically stable equilibrium point. Then as the time constant of the inhibitory postsynaptic potential increases, an effect of certain anesthetic agents on the single neuron, this equilibrium point becomes unstable, giving rise to a stable limit cycle. With further increases in the inhibitory time constant the limit cycle disappears and we again find an asymptotically stable equilibrium point corresponding to a decreasing mean excitatory synaptic drive. This offers a new interpretation of the stochastic characteristic of the induction of anesthesia. There have been animal studies, and also a small number of human studies, in which the response of a single subject to varying doses of an anesthetic agent are observed. In these studies, it has been observed that the concentrationresponse curve is very sharp [30]. In particular, as the dose or concentration of the anesthetic agent increases there exists a point at which there is an abrupt transition from consciousness to unconsciousness. However, in these studies the concentrationresponse curve has not been found to be a strict step function. Rather, there is a narrow range of concentrations within which the therapeutic effect (response or lack of a response to a noxious stimulus) has been described as a probability of response.
A common empirical model describing this effect is a variation on the Hill equation [31] in which the probability P of an anesthetic effect, that is, of no response to a noxious stimulus, is given by where C is the concentration of anesthetic agent, C 50 is the concentration of anesthetic agent associated with a 50 % probability of an anesthetic effect, and γ is a parameter that determines the steepness of the concentration-probability anesthetic curve. Typically for individual subjects γ ranges from 6 to 20 [32,33]. The range of drug concentrations where we cannot predict with certainty whether the patient will respond to noxious stimuli or not (requiring a probabilistic model) could certainly be explained by stochastic variation in the postulated inputs from pain, sensorimotor, and proprioceptive sensors, or the reticular activating system. However, an alternative explanation is that during the induction of anesthesia the probability of observing the response of the patient depends on the time in the cycle at which the observation is made. The patient will respond if the observation occurs at the time when the mean excitatory synaptic drive is at its peak and will not respond if the mean excitatory synaptic drive is at its nadir. If the concentration of the anesthetic drug is further increased, corresponding to a further increase in the inhibitory postsynaptic potential time constant, the limit cycle disappears and the mean excitatory synaptic drive has an asymptotically stable equilibrium that approaches zero (Fig. 10) and at this point, the patient is deeply anesthetized (i.e., sedated) and will not respond.
The most formidable challenge to understanding the transition to anesthesia is the immense dimensionality of the brain. Typically, this dimensionality is simplified with either specific structural assumptions or by using mean field (or more sophisticated statistical) frameworks. The extent to which these assumptions influence even qualitative conclusions is unclear. In the case of the mean synaptic drive model that we use in this paper, we have simplified the analysis by ignoring second-order terms in an expansion around mean values. Without this assumption, we are unable to reach any analytical conclusions regarding excitatory neuronal activity with an increasing inhibitory postsynaptic potential time constant. It is not clear how the mean field approximation impacts our conclusions.
In the second part of the paper we have accordingly extended our analysis without mean field assumptions by postulating the existence of a subset of inhibitory neurons that themselves do not receive inhibitory inputs. Using this assumption we derived sufficient conditions for convergence and state equipartition of the activity of excitatory neurons receiving input from these inhibitory neurons to zero. We further extended this analysis to consider a nonmonotonic postsynaptic potential, as an alternative and possibly more realistic model for postsynaptic potentials. An assumption that is essential to this finding of our paper is the existence of a subset of inhibitory neurons that themselves do not receive inhibitory input. We view this postulate as a more transparent assumption for the synaptic drive model than the mean field assumption. We are unaware of any experimental verification of such a subset. In the absence of such data, we view this assumption as defining a sufficient condition for partial state equipartitioning of a subset of excitatory neurons to zero activity. The convergence of activity of subset(s) of excitatory neurons would certainly be a plausible explanation for the interruption of information flow in the brain that is one suggested mechanism for unconsciousness [20]. The proposed dynamical system framework can potentially foster the development of new frameworks that can allow us to interpret experimental and clinical results, connect biophysical findings to psychophysical phenomena, explore a new hypothesis based on the cognitive neuroscience of consciousness and develop new assertions, and improve the reliability of general anesthesia. Dynamical systems theory is ideally suited for rigorously describing the behavior of large-scale networks of neurons and can potentially establish a scientific basis for new metrics of anesthetic depth by making the assessment of consciousness a mechanistically grounded tool.