Human Brain Networks: Spiking Neuron Models, Multistability, Synchronization, Thermodynamics, Maximum Entropy Production, and Anesthetic Cascade Mechanisms

Advances in neuroscience have been closely linked to mathematical modeling beginning with the integrate-and-fire model of Lapicque and proceeding through the modeling of the action potential by Hodgkin and Huxley to the current era. The fundamental building block of the central nervous system, the neuron, may be thought of as a dynamic element that is “excitable”, and can generate a pulse or spike whenever the electrochemical potential across the cell membrane of the neuron exceeds a threshold. A key application of nonlinear dynamical systems theory to the neurosciences is to study phenomena of the central nervous system that exhibit nearly discontinuous transitions between macroscopic states. A very challenging and clinically important problem exhibiting this phenomenon is the induction of general anesthesia. In any specific patient, the transition from consciousness to unconsciousness as the concentration of anesthetic drugs increases is very sharp, resembling a thermodynamic phase transition. This paper focuses on multistability theory for continuous and discontinuous dynamical systems having a set of multiple isolated equilibria and/or a continuum of equilibria. Multistability is the property whereby the solutions of a dynamical system can alternate between two or more mutually exclusive Lyapunov stable and convergent equilibrium states under asymptotically slowly changing inputs or system parameters. In this paper, we extend the theory of multistability to continuous, discontinuous, and stochastic nonlinear dynamical systems. In particular, Lyapunov-based tests for multistability and synchronization of dynamical systems with continuously differentiable and absolutely continuous flows are established. The results are then applied to excitatory and inhibitory biological neuronal networks to explain the underlying mechanism of action for anesthesia and consciousness from a multistable dynamical system perspective, thereby providing a theoretical foundation for general anesthesia using the network properties of the brain. Finally, we present some key emergent properties from the fields of thermodynamics and electromagnetic field theory to qualitatively explain the underlying neuronal mechanisms of action for anesthesia and consciousness.

Abstract: Advances in neuroscience have been closely linked to mathematical modeling beginning with the integrate-and-fire model of Lapicque and proceeding through the modeling of the action potential by Hodgkin and Huxley to the current era. The fundamental building block of the central nervous system, the neuron, may be thought of as a dynamic element that is "excitable", and can generate a pulse or spike whenever the electrochemical potential across the cell membrane of the neuron exceeds a threshold. A key application of nonlinear dynamical systems theory to the neurosciences is to study phenomena of the central nervous system that exhibit nearly discontinuous transitions between macroscopic states. A very challenging and clinically important problem exhibiting this phenomenon is the induction of general anesthesia. In any specific patient, the transition from consciousness to unconsciousness as the concentration of anesthetic drugs increases is very sharp, resembling a thermodynamic phase transition. This paper focuses on multistability theory for continuous and discontinuous dynamical systems having a set of multiple isolated equilibria and/or a continuum of equilibria. Multistability is the property whereby the solutions of a dynamical system can alternate between two or more mutually exclusive Lyapunov stable and convergent equilibrium states under asymptotically slowly changing inputs or system parameters. In this paper, we extend the theory of multistability to continuous, discontinuous, and stochastic nonlinear dynamical systems. In particular, Lyapunov-based tests for multistability and synchronization of dynamical systems with continuously differentiable

Introduction
Advances in neuroscience have been closely linked to mathematical modeling beginning with the integrate-and-fire model of Lapicque [1] and proceeding through the modeling of the action potential by Hodgkin and Huxley [2] to the current era of mathematical neuroscience; see [3,4] and the numerous references therein. Neuroscience has always had models to interpret experimental results from a high-level complex systems perspective; however, expressing these models with dynamic equations rather than words fosters precision, completeness, and self-consistency. Nonlinear dynamical system theory, in particular, can provide a framework for a rigorous description of the behavior of large-scale networks of neurons. A particularly interesting application of nonlinear dynamical systems theory to the neurosciences is to study phenomena of the central nervous system that exhibit nearly discontinuous transitions between macroscopic states. One such example exhibiting this phenomenon is the induction of general anesthesia [5][6][7][8][9].
The rational, safe, and effective utilization of any drug in the practice of medicine is grounded in an understanding of the pharmacodynamics of the drug, loosely defined as what the drug does to the body [10]. A very important measure of the pharmacodynamics of any drug is the drug concentration parameter EC 50 , which reflects the drug dose at which the therapeutic effect is achieved in 50% of the cases. This concept is certainly applicable for the administration of general inhalational anesthetics, where the potency of the drug is defined by the minimum alveolar concentration (MAC) of the drug needed to prevent a response to noxious stimuli in 50% of administrations [11].
The MAC concept is intrinsically embedded in a probabilistic framework [10]. It is the concentration at which the probability of a response to a noxious stimulus is 0.5. Typically the MAC of a particular anesthetic is determined by administering various doses of the agent to a population of patients and determining the dose at which there is a 0.5 chance of responding to a noxious stimulus. (Technically, we identify the concentration in the alveoli, the fundamental functional gas exchange units of the lung, at which the chance of response is 0.5.) It has been possible, however, to conduct studies of single subjects, varying the anesthetic concentration and determining responsiveness. When this has been done, it has been noted that the transition from responsiveness to non-responsiveness in the individual patient is very sharp, almost an all-or-none transition [12]. This simply confirms the observations of generations of clinicians. And this raises the question of how to account for such a transition in terms of the known molecular properties of the anesthetic agent.
The mechanism of general anesthesia is still under considerable investigation. Theories range from a nonspecific perturbation of the lipid bilayer membrane of neurons, the cells responsible for the "information" function of the central nervous system, to the interaction of the anesthetic agent with specific protein receptors [13]. It is certainly possible that if the mechanism of general anesthesia is the binding of the anesthetic agent to a specific receptor protein, then the nearly all-or-none transition or bifurcation from the awake state to the anesthetized state could be explained by a highly cooperative binding of the anesthetic to the receptor. In fact, it has been common to mathematically model the probability of responsiveness to drug concentration using the Hill equation, a simplified static equation originally derived in 1909 to describe the cooperative binding of oxygen to the hemoglobin molecule [14]. However, an alternative explanation could be sought in the dynamic network properties of the brain.
The human central nervous system involves a complex large-scale interconnected neural network involving feedforward and feedback (or recurrent) networks, with the brain serving as the central element of this network system. The brain is interconnected to receptors that transmit sensory information to the brain, and in turn the brain delivers action commands to effectors. The neural network of the brain consists of approximately 10 11 neurons (nerve cells) with each having 10 4 to 10 5 connections interconnected through subnetworks or nuclei. The nuclei in turn consist of clusters of neurons each of which performs a specific and defined function.
The most basic characteristic of the neurons that comprise the central nervous system is the electrochemical potential gradient across the cell membrane. All cells of the human body maintain an electrochemical potential gradient between the inside of the cell and the surrounding milieu. Neurons have the capacity of excitability. If stimulated beyond a threshold, then the neuron will "fire" and produce a large voltage spike (the action potential) before returning to the resting potential [3,4]. The neurons of the brain are connected in a complex network in which the firing of one neuron can be the stimulus for the firing of another neuron.
A major focus of theoretical neuroscience has been describing neuronal behavior in terms of this electrochemical potential, both at the single neuron level but more ambitiously, at the level of multi-neuron networks. In this type of analysis the specific properties of the single neuron that are most relevant are how the spike of a one neuron alters the electrochemical potential of another neuron, and how this change in the potential results in a neuronal spike. The physical connection between neurons occurs in the synapse, a small gap between the axon, the extension of the cell body of the transmitting neuron, and the dendrite, the extension of the receiving neuron. The signal is transmitted by the release of a neurotransmitter from the axon into the synapse. This neurotransmitter diffuses across the synapse, binds to a postsynaptic receptor membrane protein on the dendrite, and alters the electrochemical potential of the receiving neuron.
It is possible that the anesthetic bifurcation to unconsciousness or the nearly all-or-none characteristic induction of anesthesia is a type of phase transition of the neural network. This possibility was first considered by Steyn-Ross et al. (see [15] and the references therein). Their focus was on the mean voltage of the soma, or cell body, of neurons. Specifically, the authors in [15] show that the biological change of state to anesthetic unconsciousness is analogous to a thermodynamic phase change involving a liquid to solid phase transition. For certain ranges of anesthetic concentrations, their first-order model predicts the existence of multiple steady states for brain activity leading to a transition from normal levels of cerebral cortical activity to a quiescent, low-firing state.
In this paper, we develope an alternative approach to the possibility of neuronal network phase transition in terms of neuronal firing rates, using the concept of multistability for dynamical systems. Multistability is the property whereby the solutions of a dynamical system can alternate between two or more mutually exclusive Lyapunov stable and convergent states under asymptotically slowly changing inputs or system parameters. In particular, multistable systems give rise to the existence of multiple (isolated and/or a continuum of) Lyapunov stable equilibria involving a quasistatic-like behavior between these multiple semistable steady states [16][17][18]. Semistability is the property whereby the solutions of a dynamical system converge to Lyapunov stable equilibrium points determined by the system initial conditions [19,20].
Multistability is ubiquitous in biological systems ranging from biochemical networks to ecosystems to gene regulation and cell replication [21][22][23]. Since molecular studies suggest that one possible mechanism of action of anesthetics is the inhibition of synaptic transmission in cortical neurons [24,25], this suggests that general anesthesia is a phenomenon in which different equilibria can be attained with changing anesthetic agent concentrations. Hence, multistability theory can potentially provide a theoretical foundation for describing general anesthesia.
Although general anesthesia has been used in the clinical practice of medicine for over 150 years, the mechanism of action is still not fully understood [13] and is still under considerable investigation [5][6][7][8][9]. Early theories postulated that anesthesia is produced by disturbance of the physical properties of cell membranes. The work of Meyer and Overton [26,27] demonstrated that for some anesthetics there was a correlation between anesthetic potency and solubility in fat-like solvents. This led to a theory that anesthesia resulted from a nonspecific perturbation of the lipid bilayer membrane of neurons [9,28]. Subsequent research then found that membrane proteins performed functions of excitability and this led to a focus on anesthetic binding and perturbation of hydrophobic regions of membrane proteins [29]. Further research also revealed that some anesthetic gases follow the Meyer-Overton correlation but do not produce anesthesia and some Meyer-Overton gases are excitatory and can cause seizures [30,31]. These results led to the more common modern focus on the interaction of the anesthetic agent with specific protein receptors [13].
In particular, there has been extensive investigation of the influence of anesthetic agents on the binding of neurotransmitters to their postsynaptic receptors [7][8][9]. A plethora of receptors have been investigated, including receptors for glycine, serotonin type 2 and 3, N-methyl-d-aspartate (NMDA), α-2 adrenoreceptors, α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA), histamine, acetylcholine, and γ-aminobutyric acid (GABA). One attractive aspect of this focus on postsynaptic receptors is it facilitates mathematical analysis on the basis of the effect of receptor binding on the postsynaptic potential. This is in marked contrast to the Meyer-Overton hypothesis, which failed to explicitly detail how a nonspecific perturbation of the lipid membrane would result in the anesthetic state.
In parallel with the investigation of the molecular interactions of general anesthetic agents, there has also been active investigation of the anatomic pathways involved in the transition from consciousness to anesthesia [5]. There is compelling evidence that the immobility created by some anesthetics is mediated at the level of the spinal cord. In contrast, functional imaging and electroencephalograph analysis has suggested that the site of suppression of consciousness is the thalamus, and thalamocortical tracts may play a critical role in the suppression of consciousness [9].
Despite these advances in our understanding of the molecular interactions of anesthetic agents and of specific anatomic loci for the action of anesthetic agents, there has been less development of a mathematical framework to understand this fascinating and clinically important phenomenon. It is certainly possible that if the mechanism of general anesthesia is the binding of the anesthetic agent to a specific receptor protein, then the nearly all-or-none transition from the awake state to the anesthetized state could be explained by a highly cooperative binding of the anesthetic to the receptor. In fact, as noted above, it has been common to mathematically model the probability of responsiveness to drug concentration using the Hill equation [14]. However, to date, no single unifying receptor mediating general anesthesia has been identified.
Rather, the most likely explanation for the mechanisms of action of anesthetics lies in the network properties of the brain. It is well established that there are two general types of neurons in the central nervous system-excitatory and inhibitory-interconnected in a complex dynamic network. The action potential of a spiking neuron is propagated along the axon to synapses where chemical neurotransmitters are released that generate a postsynaptic potential on the dendrites of connected neurons. Excitatory neurons generate a depolarizing postsynaptic potential on the dendrite of the connected neuron and if the depolarization is of sufficient magnitude, then a spike will be induced in the connected neuron. In contrast, inhibitory neurons generate a hyperpolarizing postsynaptic potential; an effect that acts to maintain a quiescent state.
There is considerable evidence that general anesthetics alter postsynaptic potentials [24,25]. An interesting example of how changes in the postsynaptic potential may be applied to the analysis of the induction of anesthesia is the view of anesthesia as a phase transition proposed by Steyn-Ross et al. (see [15] and the references therein). While their analysis was highly informative, in this paper we use a dynamical system theory framework in terms of neuronal firing rates, using the concepts of network thermodynamics [32] and multistability, for explaining the mechanisms of action for general anesthesia. This facilitates a focus on the network properties of the central nervous system. The firing rate models used for network analysis must have sufficient generality and include parameters that can account for such relevant physiological changes at the single neuron level. The synaptic drive firing model, introduced by Ermentrout and his collaborators [4,32], and the system thermodynamic framework, introduced in [32], is the underlying framework for this paper.
In this paper, the fundamental building block of the central nervous system, the neuron, is represented as a dynamic element that is "excitable", and can generate a pulse or spike whenever the electrochemical potential across the cell membrane of the neuron exceeds a threshold value. More specifically, a nonlinear discontinuous system framework is developed in Sections 2 and 3 for describing the relationship between the synaptic voltage and firing rates of excitatory and inhibitory neural networks. To establish convergence and semistability for discontinuous dynamical systems we introduce the notion of nontangency between a discontinuous vector field and a weakly invariant or weakly negatively invariant subset of level or sublevel sets of Lyapunov functions in Section 5. Specifically, to capture the notion of nontangency we introduce the direction cone of a discontinuous vector field. Then, using positive limit sets, restricted prolongations, and nontangency we develop Lyapunov analysis for convergence and semistability to establish multistability for discontinuous dynamical systems. Here, the restricted prolongation of a point is a subset of its positive prolongation as defined in [33]. In addition, using nontangency, we present Lyapunov results for convergence and semistability to develop sufficient conditions for multistability for discontinuous dynamical systems.
While previous treatments of nontangency-based Lyapunov tests for convergence and semistability for dynamical systems with continuous vector fields are given in [19], our results involve dynamical systems with discontinuous vector fields for capturing plasticity (i.e., dynamic network connections) in our neural network model generating absolutely continuous solutions necessitating stability analysis via nonsmooth Lyapunov stability theory involving Clarke generalized gradients and set-valued Lie derivatives. Using the aforementioned dynamical system framework, we apply the results of Sections 4 and 5 to excitatory-inhibitory firing neural models in an attempt to understand the mechanisms of action of anesthetics. While there is ongoing debate as to whether information is encoded by the firing rates (i.e., rate coding) of spiking neurons or by precise timing of single neuron spikes (i.e., temporal coding) [34], it is evident that firing rates do characterize central nervous system activity. Firing rates are nonnegative entities and the nonnegativity constraint for neural network activity can be easily incorporated within nonlinear dynamical system theory using solutions of differential equations and differential inclusions evolving in cones [35].
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. This may be particularly relevant to understand the mechanism of action for general anesthesia. In Section 7, we provide sufficient conditions for global asymptotic and exponential synchronization for our excitatory and inhibitory cortical neuronal network.
To avoid the complexity of large-scale and high connectivity models of the neural network of the brain, the scale and connectivity of the network can be simplified using mean field theories. Early mean field theories assumed that the brain is organized into a limited number of pools of identical spiking neurons [36]. However, more commonly mean field theories assume that the strength of connection between neurons is normally distributed around some mean value. Mean field theories can impose self-consistency on field variables; for example, if postsynaptic potentials are assumed to be a function of some mean firing rate, then those postsynaptic potentials should lead to a consistent predicted mean firing rate. The idea of applying mean field theories, drawn from the study of condensed matter, originated with [37]. Subsequently, Sompolinsky et al. [38] developed a mean field theory for neural networks analogous to the equations developed for spin glasses with randomly symmetric bonds [39]. The authors in [40] investigated the stability of system states for a network of integrate-and-fire neurons, whereas the authors in [41] extended this theoretical model to the analysis of oscillations. Gerstner et al. [42,43] subsequently developed a mean field theory using a spike response model and demonstrated that the integrate-and-fire model was a special case of the spike response model.
In Section 8 of the paper, we extend our results further by demonstrating multistability in the mean when the coefficients of the neuronal connectivity matrix are random variables. Specifically, we use a stochastic multiplicative uncertainty model to include modeling of a priori uncertainty in the coefficients of the neuronal connectivity matrix by means of state-dependent noise. Our stochastic multiplicative uncertainty model uses state-dependent Gaussian white noise to represent parameter uncertainty by defining a measure of ignorance, in terms of an information-theoretic entropy, and then determining the probability distribution which maximizes this measure subjected to agreement with a given model. To account for time delay and memory effects in inhibitory and excitatory networks, in Sections 9 and 10 we extend the results of Section 8 to a large-scale excitatory and inhibitory synaptic drive firing rate model with time-varying delays and stochastic input uncertainty, and global mean-square synchronization of this model is investigated.
Finally, in Sections 11 and 12 we discuss key emergent properties from the fields of thermodynamics and electromagnetic field theory for developing plausible mechanisms of action for general anesthesia. Specifically, in Section 11, we highlight how the supreme law of nature-the second law of thermodynamics-can be used to arrive at mechanistic models for the anesthetic cascade using the principle of maximum entropy production and thermodynamics as applied to the human brain. In Section 12, we use the idea of anesthetics disrupting the inflow of free energy to the brain as electrical signals that generate electromagnetic fields causing a shielding effect leading to the emergence of unconsciousness.

Biological Neural Networks: A Dynamical Systems Approach
The fundamental building block of the central nervous system, the neuron, can be divided into three functionally distinct parts, namely, the dendrites, soma (or cell body), and axon (see Figure 1). The dendrites play the role of input devices that collect signals from other neurons and transmit them to the soma; whereas the soma generates a signal that is transmitted to other neurons by the axon. The axons of other neurons connect to the dendrites and soma surfaces by means of connectors called synapses. The behavior of the neuron is best described in terms of the electrochemical potential gradient across the cell membrane. If the voltage gradient across the membrane increases to a critical threshold value, then there is a subsequent abrupt step-like increase in the potential gradient, the action potential. This action potential is transmitted from the soma along the axon to a dendrite of a receiving neuron. The action potential elicits the release of neurotransmitter molecules that diffuse to the dendrite of a "receiving" neuron. This alters the voltage gradient across the receiving neuron.
The electrochemical potential for a neuron can be described by a nonlinear four-state system [32]. Coupling these system equations for each neuron in a large neural population is computationally prohibitive. To simplify the mathematical modeling, it has been common to use phenomenological firing rate models for studying neural coding, memory, and network dynamics [4]. Firing rate models involve the averaged behavior of the spiking rates of groups of neurons rather than tracking the spike rate of each individual neuron cell. In such population models, the activity of a neuron, that is, the rate at which the neuron generates an action potential (i.e., "fires") is modeled as a function of the voltage (across the membrane). The "firing" of a neuron evokes voltage changes, postsynaptic potentials on receiving neurons; that is, neurons electrically connected to the firing neurons via axon-dendrite connections. means of connectors called synapses. The behavior of the neuron is best described in terms of the electrochemical potential gradient across the cell membrane. If the voltage gradient across the membrane increases to a critical threshold value, then there is a subsequent abrupt step-like increase in the potential gradient, the action potential. This action potential is transmitted from the soma along the axon to a dendrite of a receiving neuron. The action potential elicits the release of neurotransmitter molecules that diffuse to the dendrite of a "receiving" neuron. This alters the voltage gradient across the receiving neuron. The electrochemical potential for a neuron can be described by a nonlinear four-state system [36]. Coupling these system equations for each neuron in a large neural population is computationally prohibitive. To simplify the mathematical modeling, it has been common to use phenomenological firing rate models for studying neural coding, memory, and network dynamics [36]. Firing rate models involve the averaged behavior of the spiking rates of groups of neurons In general, neurons are either excitatory or inhibitory depending on whether the postsynaptic potential increases or decreases the potential of the receiving neuron. In particular, excitatory neurotransmitters depolarize postsynaptic membranes by increasing membrane potentials and can collectively generate an action potential. Inhibitory neurotransmitters hyperpolarize the postsynaptic membrane by decreasing membrane potentials, thereby nullifying the actions of excitatory neurotransmitters and in certain cases prevent the generation of action potentials.
Biological neural network models predict a voltage in the receiving or postsynaptic neuron given by where v X i (·), X ∈ {E, I}, i = 1, . . . , n E + n 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 jth 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 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 the neuronal connectivity matrix A XY , with units of volts, contains entries A XY ij , X, Y ∈ {E, I}, representing the coupling strength of the jth neuron on the ith neuron such that A XE ij > 0 and A XI ij < 0, X ∈ {E, I}, if the jth neuron is connected (i.e., contributes a postsynaptic potential) to the ith neuron, and A XY ij = 0, otherwise. Furthermore, v E thi (·) and v I thi (·) are continuous input voltages characterizing nerve impulses from sensory (pain) receptors, sensorimotor (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 by and assuming an exponential decay of the synaptic voltages of the form ( [4,44]) where the dimensionless gain B (E,I) i is equal to B E i if the ith neuron is excitatory and B I i if the ith neuron is inhibitory, and similarly for S , it follows from Equations (4) and (5) that Now, using the expressions for the excitatory and inhibitory voltage given by Equations (2) and (3), respectively, it follows that The above analysis reveals that a form for capturing the neuroelectic 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 voltage to the ith neuron, A ij is a constant representing the coupling strength of the jth neuron on the ith neuron, τ i is a time constant, 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.
In this paper, we will explore different activation functions including discontinuous hard-limiter activation functions, and continuous half-wave rectification, saturation, and sigmoidal functions. Specifically, for a typical neuron ( [3]) where i ∈ {1, . . . , n} and [x] + = x if x ≥ 0, and [x] + = 0, otherwise. Alternatively, we can approximate f i (x) by the smooth (i.e., infinitely differentiable) half-wave rectification function where i ∈ {1, . . . , n} and γ 0. Note that f i (x) ≈ 1 for x > 0 and f i (x) ≈ 0, x = 0. In addition, note that Equations (10) and (11) reflect the fact that as the voltage increases across the membrane of the ith neuron, the firing rate increases as well. Often, the membrane potential-firing rate curve exhibits a linear characteristic for a given range of 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 where i ∈ {1, . . . , n}, γ 0, and f max = lim γ→∞ f i (x) denotes the maximum firing rate.

Connections to Mean Field Excitatory and Inhibitory Synaptic Drive Models
The excitatory and inhibitory neural network model given by Equations (7) and (8) can possess multiple equilibria. For certain values of the model parameters it can be shown that as the inhibitory time constants λ I i , i = 1, . . . , n I , get larger, multiple stable and unstable state equilibria can appear [45]. Since molecular studies suggest that one possible mechanism of action of anesthetics is the prolongation of the time constants of inhibitory neurons [24,25], this suggests that general anesthesia is a phenomenon in which different equilibria can be attained with changing anesthetic agent concentrations. To explore this multistability phenomenon, in [45,46] we developed a simplified scale and connectivity neural network model using a mean field theory. As noted in the Introduction, mean field theories assume that the brain is organized into limited number of pools of identical spiking neurons [36]. However, more commonly, mean field theories assume that the strength of the connection between neurons is normally distributed around some mean value.
To see how our general excitatory and inhibitory synaptic drive model given by Equations (7) and (8) can be reduced to a mean excitatory and mean inhibitory model, consider Equations (7) and (8) with In this case, Equations (7) and (8) become where f (·) is given by either Equation (11) or Equation (12) and ij , X, Y ∈ {E, I}, denote mean and ∆ XY ij , X, Y ∈ {E, I}, are deviations from the mean. In this case, it follows that Using the average and perturbed expression for A XY ij , X, Y ∈ {E, I}, Equations (13) and (14) can be rewritten as where S denote the mean excitatory synaptic drive and mean inhibitory synaptic drive in dimensionless units, respectively. Now, defining where δ E i (t) and δ I i (t) are deviations from the mean, Equations (16) and (17) 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 Equations (18) and (19) are small relative to the remaining terms in f (·). Then, a first-order expansion of Equations (18) and (19) gives Now, assuming that the higher-order terms can be ignored, Equations (20) and (21) become Finally, summing Equations (22) and (23) over i = 1, . . . , n E and i = 1, . . . , n I , dividing by n E and n I , respectively, using Equation (15), and assuming v E it follows that the average excitatory synaptic drive and the average inhibitory synaptic drive are given by dS Equations (24) and (25) represent the spatial average (mean) dynamics of the system given by Equations (13) and (14), and 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 [36]. Preliminary results using the simplified two-state mean excitatory and inhibitory synaptic drive model given by Equations (24) and (25) for connecting notions of general anesthesia to multistability and bifurcations are given in [45,46].

Multistability Theory and Discontinuous Spiking Neuron Models
Multistability is the property whereby the solutions of a dynamical system can alternate between two or more mutually exclusive semistable states under asymptotically slowly changing inputs or system parameters. In particular, the state of a multistable system converges to Lyapunov stable equilibria that belong to an equilibrium set that has a multivalued hybrid topological structure consisting of isolated points and closed sets homeomorphic to intervals on the real line. To develop the notion of multistability, consider the autonomous differential equation given bẏ where, for every t ≥ 0, x(t) ∈ D ⊆ R n , f : R n → R n is Lebesgue measurable and locally essentially bounded [47], that is, f is bounded on a bounded neighborhood of every point x, excluding sets of measure zero. Furthermore, let E e {x ∈ R n : f (x) = 0} denote the set of equilibria for Equation (26).
where the Filippov set-valued map K[f ] : R n → 2 R n is defined by and where B δ (x) denotes the open ball centered at x with radius δ, 2 R n denotes the collection of all subsets of R n , µ(·) denotes the Lebesgue measure in R n , "co" denotes the convex closure, and µ(S)=0 denotes the intersection over all sets S of Lebesgue measure zero [49].
Note that since f is locally essentially bounded, K[f ](·) is upper semicontinuous and has nonempty, compact, and convex values. Thus, Filippov solutions are limits of solutions to Equation (26) with f averaged over progressively smaller neighborhoods around the solution point, and hence, allow solutions to be defined at points where f itself is not defined. Hence, the tangent vector to a Filippov solution, when it exists, lies in the convex closure of the limiting values of the system vector field f (·) in progressively smaller neighborhoods around the solution point. Note that K[f ] : R n → 2 R n is a map that assigns sets to points.
Dynamical systems of the form given by Equation (27) are called differential inclusions [50] and for each state x ∈ R n , they specify a set of possible evolutions rather than a single one. It follows from 1) of Theorem 1 of [51] that there exists a set N f ⊂ R n of measure zero such that, for every set W ⊂ R n of measure zero, where {x i } i∈Z + ⊂ R n converges to x ∈ R n . Differential inclusions include piecewise continuous dynamical systems as well as switched dynamical systems as special cases. For example, if f (·) is piecewise continuous, then Equation (26) can be represented as a differential inclusion involving the Filippov set-valued map of piecewise continuous vector fields given by where S f has measure zero and denotes the set of points where f is discontinuous [52]. Similarly, differential inclusions can include, as a special case, switched dynamical systems of the forṁ where x(t) ∈ D ⊆ R n , f p : R n → R n is locally Lipschitz continuous, and p ∈ P = {1, . . . , d} is a finite index set. To see how the state-dependent switched dynamical system given by Equation (29) can be represented as a differential inclusion involving Filippov set-valued maps, define the switched system (29) with a piecewise linear partitioned state space as the triple (D, Q, V), where D = R n or D is a polytope in R n of dimension dim(D) = n, Q = {D p } p∈P is a piecewise linear partition of D with index set P, and V = {f p } p∈P n , and where f p : U p → R n , U p is an open neighborhood of D p , and P n = {p ∈ P : dim(D p ) = n}. Specifically, ∪ d p=1 D p = D ⊆ R n , where D is a polytope (i.e., the convex hull of finitely many points, and hence, compact) and D p , p = 1, . . . , d, is a family of polyhedral sets in R n with nonempty interiors. Furthermore, for every i, j ∈ {1, . . . , d}, i = j, let D i ∩ D j = ∅ or D i ∩ D j is a (n − 1)-dimensional manifold included in the boundaries ∂D i and ∂D j . Finally, since each vector field f i is Lipschitz continuous in the state x, it defines a continuously differentiable flow ψ i (t, x) within every open set U i ⊃ D i . In particular, each flow ψ i (t, x) is well defined on both sides of the boundary ∂D j . Thus, in the interior of each operating region D p , the global dynamics of Equation (29) is completely described by the local dynamics characterized by a particular vector field f p , and hence, there exists a unique classical (i.e., continuously differentiable) solution to Equation (29). However, for a point x ∈ H p ∩ D p , where H p is some supporting hyperplane of D p , nonuniqueness and nonexistence of solutions to Equation (29) can occur.
To address the problem of existence and extendability of solutions for Equation (29), let the global dynamics of Equation (29) be characterized by one of the following differential inclusionṡ Note that F : R n → 2 R n is nonempty and finite, and hence, compact. Moreover, since each f p : is convex, and hence, F c : R n → 2 R n is an upper semicontinuous set-valued map with nonempty, convex, and compact values. That is, for every x ∈ D and every ε > 0, there exists δ > 0 such that, for all z ∈ R n satisfing z − x ≤ δ, F c (z) ⊆ F c (x) + εB 1 (0). In this case, it can be shown that there exists a Filippov solution to the differential inclusion given by Equation (30) at all interior points x ∈ D ⊆ R n [47,53]. In addition, if for each unbounded cell where T x D denotes the tangent cone to D at x (see Definition 5) [54], then there exists a Filippov solution x : [0, ∞) → R n to Equation (30) for every x 0 ∈ D ⊆ R n ( [50], p. 180).
Switched dynamical systems are essential in modeling the plasticity of the central nervous system. In particular, recent neuroimaging findings have shown that the loss of top-down (feedback) processing in association with anesthetic-induced unconsciousness observed in electroencephalographic (EEG) and functional magnetic resonance imaging (fMRI) is associated with functional disconnections between anterior and posterior brain structures [55][56][57]. These studies show that topological rather than network connection strength of functional networks correlate with states of consciousness. Hence, changes in network topology are essential in capturing the mechanism of action for consciousness rather than network connectivity strength during general anesthesia. Dynamic neural network topologies capturing link separation or creation, which can be modeled by switched systems and differential inclusions, can provide a mechanism for the objective selective inhibition of feedback connectivity in association with anesthetic-induced unconsciousness.
Since We say that a set M is weakly positively invariant (resp., strongly positively invariant) with respect to Equation (26) if, for every x 0 ∈ M, M contains a right maximal solution (resp., all right maximal solutions) of Equation (26) [48,59]. The set M ⊆ R n is weakly negatively invariant if, for every x ∈ M and t ≥ 0, there exist z ∈ M and a Filippov solution ψ(·) to Equation (26) The set M ⊆ R n is weakly invariant if M is weakly positively invariant as well as weakly negatively invariant. Finally, an equilibrium point of Equation (26) is a point x e ∈ R n such that 0 ∈ K[f ](x e ). It is easy to see that x e is an equilibrium point of Equation (26) if and only if the constant function x(·) = x e is a Filippov solution of Equation (26). We denote the set of equilibrium points of Equation (26) To develop Lyapunov theory for nonsmooth dynamical systems of the form given by Equation (26), we need the notion of generalized derivatives and gradients. In this paper, we focus on Clarke generalized derivatives and gradients [52,60].
The Clarke generalized gradient ∂V : where co denotes the convex hull, ∇ denotes the nabla operator, N is the set of measure zero of points where ∇V does not exist, S is any subset of R n of measure zero, and the increasing unbounded sequence Note that Equation (32) always exists. Furthermore, note that it follows from Definition 2 that the generalized gradient of V at x consists of all convex combinations of all the possible limits of the gradient at neighboring points where V is differentiable. In addition, note that since V (·) is Lipschitz continuous, it follows from Rademacher's theorem ( [61](Theorem 6, p. 281)) that the gradient ∇V (·) of V (·) exists almost everywhere, and hence, ∇V (·) is bounded. Specifically, for every x ∈ R n , every ε > 0, and every Lipschitz constant L for V on B ε (x), ∂V (x) ⊆ B L (0). Thus, since for every x ∈ R n , ∂V (x) is convex, closed, and bounded, it follows that ∂V (x) is compact.
In order to state the results of this paper, we need some additional notation and definitions. Given a locally Lipschitz continuous function V : R n → R, the set-valued Lie derivative L f V : R n → 2 R of V with respect to f at x [48,62] is defined as If to denote the largest (resp., smallest) element of L f V (x). Furthermore, we adopt the convention The next definition introduces the notion of semistability for discontinuous dynamical systems. This stability notion is necessary for systems having a continuum of equilibria. Specifically, since every neighborhood of a nonisolated equilibrium contains another equilibrium, a nonisolated equilibrium cannot be asymptotically stable. Hence, asymptotic stability is not the appropriate notion of stability for systems having a continuum of equilibria. Two notions that are of particular relevant to such systems are convergence and semistability. Convergence is the property whereby every system solution converges to a limit point that may depend on the system initial condition (i.e., initial anesthetic concentrations). Semistability is the additional requirement that all solutions converge to limit points that are Lyapunov stable. Semistability for an equilibrium thus implies Lyapunov stability, and is implied by asymptotic stability. Thus, semistability guarantees that small perturbations from the limiting state of unconsciousness will lead to only small transient excursions from that state of unconsciousness. It is important to note that semistability is not merely equivalent to asymptotic stability of the set of equilibria. Indeed, it is possible for a trajectory to converge to the set of equilibria without converging to any one equilibrium point as examples in [19] show. For further details, see [20].  (26) is Lyapunov stable if, for every ε > 0, there exists δ = δ(ε) > 0 such that, for every initial condition x 0 ∈ B δ (z) and every Filippov solution x(t) with the initial condition (26) is semistable if z is Lyapunov stable and there exists an open subset D 0 of D containing z such that, for all initial conditions in D 0 , the Filippov solutions of Equation (26) converge to a Lyapunov stable equilibrium point. The system given by Equation (26) is semistable with respect to D if every Filippov solution with initial condition in D converges to a Lyapunov stable equilibrium. Finally, the system given by Equation (26) is said to be globally semistable if the system given by Equation (26) is semistable with respect to R n .
Finally, we introduce the definition of multistability of the dynamical system given by Equation (26).
Definition 4 Consider the nonlinear dynamical system given by Equation (26). We say that the dynamical system given by Equation (26) is multistable if (i) there exists more than one equilibrium point of Equation (26) in R n ; (ii) all Filippov solutions to Equation (26) converge to one of these equilibrium points; and (iii) almost all Filippov solutions to Equations (26) converge to Lyapunov stable equilibria; that is, the set of initial conditions driving the Filippov solutions of Equation (26) to unstable equilibria has Lebesgue measure zero.
It is important to note that our definition of multistability is different from the definition given in [18]. Specifically, pertaining to condition (iii), the definition of multistability given in [18] requires that almost all Filippov solutions to Equation (26) converge to asymptotically stable equilibria. This key difference between the two definitions allows for the dynamical system given by Equation (26) to possess a continuum of equilibria, rather than merely isolated equilibria. As we see later, if f i is of the form given by Equation (10), then Equation (9) has a continuum of equilibria under certain conditions, and hence, Equation (26) is semistable in the sense of (iii) [63]. Hence, in this case, it is appropriate to use Definition 4 to characterize multistability.
Almost all of the existing results on multistability theory rely on linearization techniques based on the Hartman-Grobman theorem [64,65] involving the fact that the linearized system has the same topological property as the original system around a hyperbolic fixed point. When the system fixed point is not hyperbolic, however, these techniques fail to predict multistability. In this case, checking multistability becomes a daunting task. Rather than checking the transversality condition for hyperbolicity, we propose a new approach for guaranteeing multistability using equilibria-independent, semidefinite Lyapunov function methods. In particular, using the geometric structure of the vector field f for a given dynamical system, we develop nontangency-based Lyapunov tests for verifying conditions (ii) and (iii) in Definition 4 involving convergence and Lyapunov stability almost everywhere.

Direction Cones, Nontangency, Restricted Prolongations, and Nonsmooth Multistability Theory
To develop multistability theory for discontinuous dynamical systems of the form given by Equation (26) we use the notions of direction cones, nontangency, and restricted prolongations. In particular, to show condition (ii) in Definition 4 holds for dynamical systems of the form given by Equation (26), we adopt the notion of nontangency [19,63] to develop nontangency-based Lyapunov tests for convergence. Specifically, the authors in [19] develop a general framework for nontangency-based Lyapunov tests for the convergence of dynamical systems described by ordinary differential equations with continuous vector fields. In [63], the authors extend some of the results of [19] to nonsmooth dynamical systems, that is, systems described by ordinary differential equations with the discontinuous right-hand sides. Since the vector field f characterizing biological neural networks can involve either continuous (e.g., half-wave rectification functions) or discontinuous (e.g., hard-limiter activation functions) vector fields, and, more importantly, the fact that anesthetics reconfigure the topological structure of functional brain networks [55][56][57] by suppressing midbrain/pontine areas involved with regulatory arousal leading to a dynamic network topology, we use the more general definition for nontangency presented in [63].
Intuitively, a vector field is nontangent to a set at a point if the vector field at the point is not contained in the tangent space to the set at that point. We use this intuitive idea for the case where the vector field describing the system dynamics is discontinuous and the set is the set of equilibria of the system. However, this notion presents two key difficulties when the set under consideration is the set of singular points of the vector field, that is, the set of equilibria of the system. In particular, the vector field at an equilibrium point is zero, and hence, it is always contained in the tangent space to the set of equilibria. In this case, in order to capture the notion of nontangency, we introduce the direction cone of a vector field.
Alternatively, the set of equilibria may not be sufficiently regular to possess a tangent space at the equilibrium point under consideration and may have corners or self-intersections. For example, in firing rate population models appearing in neuroscience, the firing rate is a nonnegative quantity representing the probability of the firing action potential by the neuron and can be interpreted as a measure of the neuron's activity. Since the firing rate of the excitatory-inhibitory network is nonnegative, all solutions of physical interest always take values in the nonnegative orthant R n + of the state space for nonnegative initial conditions. For such systems, which evolve on possibly closed positively invariant subsets of R n , it is natural to consider the nonnegative orthant as their state space. Hence, the dynamical system evolves on the nonnegative orthant and can have the boundary of the orthant as its set of equilibria. In this case, the set of equilibria has a corner at the origin. We overcome this difficulty by considering the tangent cone [66,67], which extends the notion of a tangent space to a nonsmooth setting.
To introduce the notions of direction cone and tangent cone, some notation and definitions are required. A set E ⊆ R n is connected if and only if every pair of open sets U i ⊆ R n , i = 1, 2, satisfying E ⊆ U 1 ∪ U 2 and U i ∩ E = ∅, i = 1, 2, has a nonempty intersection. A connected component of the set E ⊆ R n is a connected subset of E that is not properly contained in any connected subset of E. Given a set E ⊆ R n , let coco E denote the convex cone generated by E.
The notion of nontangency introduced in Definition 5 is different from the well-known notion of transversality [68]. Transversality between a vector field and a set is possible only at a point in the set where the vector field is not zero and the set is locally a differentiable submanifold of codimension one. Alternatively, nontangency is possible even if the vector field is zero and the set is not a differentiable submanifold of codimension one.
Definition 5 formalizes the notion of nontangency by defining nontangency of a discontinuous vector field to a set at a point to be the condition that the tangent cone to the set at the point and the direction cone of the vector field at that point have no nonzero vector in common. Using the notion of nontangency, in [63] we developed necessary and sufficient conditions for convergence of discontinuous dynamical systems. Specifically, convergence of system trajectories are guaranteed if and only if the vector field is nontangent to the positive limit set of the point at some positive limit point. However, this result cannot be applied directly in practice since it is not generally possible to find the positive limit set of system solution trajectories. Since nontangency to any outer estimate of the positive limit set implies nontangency to the positive limit set itself, we use nontangency-based Lyapunov tests for convergence. In particular, if the vector field f is nontangent to the largest invariant subset of the zero-level set of the derivative of a Lyapunov function that is nonincreasing along the solutions of Equation (26), then every bounded solution converges to a limit.
Since the application of the convergence results discussed above depends on verifying the boundedness of trajectories, the well-known results for boundedness involving proper (that is, radially unbounded in the case where the state space is R n ) Lyapunov functions [69,70] by introducing the notion of a weakly proper function need to be extended. Specifically, in [63] we consider Lyapunov functions whose connected components of their sublevel sets are compact. In this case, the existence of a weakly proper Lyapunov function that is nonincreasing along the system trajectories implies that the trajectories are bounded.
Using the notion of nontangency we then developed Lyapunov measures for almost everywhere semistability to arrive at multistability theory for discontinuous dynamical systems [63]. Here, prolongations [33,71] play a role analogous to that played by positive limit sets in the aforementioned discussion. In particular, the notion of a restricted prolongation of a point is used to show that an equilibrium point of Equation (26) is Lyapunov stable if and only if the discontinuous vector field is nontangent at the equilibrium to its restricted prolongation. The restricted prolongation of a point is a subset of its positive prolongation [33,71] and is defined as follows.
Definition 6 Given a point x ∈ R n and a bounded open neighborhood U ⊂ R n of x, the restricted prolongation of x with respect to U is the set R U x ⊆ U of all subsequential limits of sequences of the form As in the case for positive limit sets discussed above, since it is not generally possible to find the restricted prolongation of an equilibrium point in practice and since nontangency to any outer estimate of the restricted prolongation implies nontangency to the restricted prolongation itself, we use outer estimates of restricted prolongations in terms of connected components of invariant and negatively invariant subsets of level and sublevel sets of Lyapunov functions and their derivatives. By assuming nontangency of the vector field to invariant or negatively invariant subsets of the level set of the Lyapunov function containing the equilibrium we can trap the restricted prolongation and the positive limit set, respectively, in the level sets of the Lyapunov function and its derivative.
The following theorem establishes sufficient conditions for convergence and Lyapunov stability almost everywhere for the system given by Equation (26). This result follows from the fact that R U x is connected and the fact that if N is composed of isolated equilibria of Equation (26) or f is nontangent to N at every point in N , then the solution to Equation (26) is convergent. For details of these facts; see [45]. For the statement of the next result,V denotes the set-valued Lie derivative for the Filippov solutions to Equation (26) and

Theorem 1 ([45])
Assume there exists a locally Lipschitz continuous and regular function V : R n → R such thatV is defined almost everywhere on R n and satisfiesV (x) ≤ 0 for almost all x ∈ R n . Let x ∈ R n be such that the solution of Equation (26) is bounded and let N denote the largest weakly invariant set contained inV −1 (0). If either every point in N is Lyapunov stable or f is nontangent to N at every point in N , then almost all solutions of Equation (26) converge to Lyapunov stable equilibria.
Example 1 Consider the two-class mean excitatory and mean inhibitory synaptic drive network characterized by a discontinuous vector field for modeling plasticity in the network given bẏ where S E 0 ≥ 0, S I 0 ≥ 0, and f i (·), i = 1, 2, are given by where step(y) = 1 for y ≥ 0 and step(y) = 0, otherwise. For this system, the set of equilibria in Next, we show that for almost all the initial conditions (S E 0 , S I 0 ) ∈ R 2 + , the equilibrium set E is attractive. To see this, consider the function V : It can be verified thatV (S E , S I ) ≤ 0 for almost all (S E , S I ) ∈ R 2 + andV −1 (0) = E. Next, we show that the vector field f of the system given by Equations (35) and (36) is nontangent to E. Let (S E , S I ) ∈ E and note that it follows from the expression of f that the direction cone F (S E ,S I ) of the vector field f at (S E , S I ) ∈ E is given by In addition, the tangent cone to E at (S E , S I ) ∈ E is given by Now, it follows from Equations (39) and (40) that T (S E ,S I ) E ∩ F (S E ,S I ) = {0}, and hence, for every (S E , S I ) ∈ E, f is nontangent to E at (S E , S I ).
Finally, it follows from Theorem 5.1 that almost all solutions of the system converge to Lyapunov stable equilibria in E, and hence, by definition, the system given by Equations (35) and (36) is multistable. Figure 2 shows the state trajectories versus time for the initial condition [S E 0 , S I 0 ] = [1,2]. This system is additionally synchronized; a notion that we discuss in Section 7.

Multistability of Excitatory-Inhibitory Biological Networks
Having developed multistability theory for discontinuous dynamical systems, in this section we apply the results of Section 5 to the excitatory-inhibitory neural firing rate model given by Equation (9). The form of biological neural network models given by Equation (9) represents a wide range of firing rate population models appearing in neuroscience [3,4]. The firing rate is a nonnegative quantity representing the probability of the firing action potential by the neuron and can be interpreted as a measure of the neuron's activity. Since the firing rate of the excitatory-inhibitory network is nonnegative, all solutions of physical interest always take values in the nonnegative orthant of the state space for nonnegative initial conditions. For such systems, which evolve on possibly closed positively invariant subsets of R n , it is natural to consider the nonnegative orthant R n + as their state space, and hence, these systems are nonnegative dynamical systems [35]. In this case, the stability and convergence result developed in Section 5 holds with respect to R n + by replacing R n with R n + . For related details, see [35]. The following result, which follows from Proposition 2.1 of [35], gives necessary and sufficient conditions for the firing rates S i (t), i = 1, . . . , n, to remain in the nonnegative orthant of the state space. For the statement of the next result recall that f is nonnegative [35] if and only if f (x) ≥≥ 0, x ∈ R n + , where "≥≥" denotes a component-wise inequality.
Proposition 1 Consider the excitatory-inhibitory network given by Equation (9). The firing rate vector S(t) [S 1 (t), . . . , S n (t)] T ∈ R n remains in the nonnegative orthant of the state space R n + for all t ≥ 0 if and only if, for every S i ≥ 0 and v thi ≥ 0, i = 1, . . . , n, the functionf = [f 1 ( n j=1 A 1j S j + v th1 ), . . . , f n ( n j=1 A nj S j + v thn )] T : R n → R n is nonnegative.
Finally, let W be the largest weakly invariant set contained inV −1 (0) and note that since S(t) ∈ R n + , t ≥ 0, W ⊆ N (Ω 1 )∩R n + . Now, since S ∈ N (Ω 1 )∩R n + is Lyapunov stable with respect to Equation (41), it follows from Corollary 4.2 of [45] that all the solutions of the excitatory-inhibitory network given by Equation (41) converge to one of the Lyapunov stable equilibria in N (Ω 1 ) ∩ R n + for Equation (41) with S 0 ∈ R n + . Hence, it follows from Theorem 1 that Equation (41) is multistable.
Example 2 Consider the excitatory-inhibitory network characterized by the dynamicṡ where S 01 ≥ 0, S 02 ≥ 0, and, f i (·), i = 1, 2, are defined by Equation (10), and note that Equations (48) and (49)  Next, it follows from Proposition 1 that S(t) ≥≥ 0, t ≥ 0, for all (S 10 , S 20 ) ∈ R 2 + , and hence, Equations (48) and (49) collapse to the linear model given bẏ Clearly, the system given by Equations (50) and (51) is Lyapunov stable, and hence, every S ∈ N (Ω 1 ) ∩ R n + is Lyapunov stable with respect to Equations (48) and (49). Now, it follows from Theorem 2 that the system given by Equations (48) and (49) is mutistable. For the initial conditions S 01 = 3 and S 02 = 1 the trajectories of the state variables with respect to time is shown in Figure 3.

Synchronization of Biological Neural Networks
Numerous complex large-scale dynamical networks often demonstrate a degree of synchronization. System synchronization typically involves coordination of events that allows a dynamical system to operate in unison resulting in system self-organization. The onset of synchronization in populations of coupled dynamical networks have been studied for various complex networks including network models for mathematical biology, statistical physics, kinetic theory, bifurcation theory, as well as plasma physics [72]. Synchronization of firing neural oscillator populations using probabilistic analysis has also been addressed in the neuroscience literature [73]. One of the most important questions in neuroscience is how do neurons, or collections of neurons, communicate. In other words, what is the neural code? There is extensive experimental verification that collections of neurons may function as oscillators [74][75][76] and the synchronization of oscillators may play a key role in the transmission of information within the central nervous system. This may be particularly relevant to understanding the mechanism of action for general anesthesia [45].
It has been known for a long time that general anesthesia has profound effects on the spectrum of oscillations in the electroencephalograph [77,78]. More recently, the authors in [79] 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 [31].
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. To develop global synchronization properties for the biological neural network system given by Equation (41) we introduce the notions of asymptotic synchronization and exponential synchronization.

Definition 7
The biological neural network given by Equation (41) is said to be globally asymptotically synchronized if lim for all S 0 ∈ R n + and i, j = 1, 2, . . . , n, i = j.
The following theorems provide sufficient conditions for global asymptotic synchronization and global exponentially synchronization of the biological neural network system given by Equation (41). For the statement of the theorems we define the ones vector of order n by e n [1, . . . , 1] T .
Proof. The proof is similar to the proof of Theorem 3 using the function V : [0, ∞) × R n + → R given by V (t, S) = e 2εt S T P S and, hence, is omitted.

Definition 9 ([80])
Let (R, +, ·) be a ring with the two binary operations of addition (+) and multiplication (·) connected by distributive laws, and let T (R, K) be the set of matrices with entries in R such that the sum of the entries in each row is equal to K for some K ∈ R.

Lemma 1 ([81])
Let G be an n × n matrix such that G ∈ T (R, K). Then there exists a matrix G ∈ R (n−1)×(n−1) given by G = M GJ such that M G = GM , where M and J are given by Next, we analyze the synchronization properties of a special class of biological neural network given by Equation (41) where L, A ∈ T (R, K) and f i (·), i = 1, . . . , n, satisfies where f max denotes the maximum firing rate.
Lemma 2 Consider the biological neural network given by Equation (41) and assume that f i (·), i = 1, . . . , n, is given by Equation (61). Then S(t) is bounded for all t ≥ 0.
Proof. Consider the function V (S) = S T S and note that the derivative of V (S) along the trajectories of Equation (41) satisfiesV Note that since 2x T y ≤ rx T x + 1 r y T y for all x, y ∈ R n and r > 0, it follows thaṫ Now, since 0 ≤ f i (x) ≤ f max for all x ∈ R and S T LS ≥ λ min (L)S T S, it follows thaṫ Next, we show that if S 0 2 = c < m, where m = n λ 2 min (L) f 2 max , then V (S(t)) ≤ m for all t ≥ 0. To see this, assume, ad absurdum, that V (S(t)) > m for some t > 0, which holds since V (S(t)) is continuous in t, and note that there exists τ ∈ (0, t) such that V (S(τ )) = m andV (S(τ )) > 0. Now, using Equation (64) it follows thatV (S(τ )) ≤ 0 for all τ > 0 such that V (S(τ )) = m, which contradictsV (S(τ )) > 0. Alternately, if S 0 2 = c ≥ m, then using a similar argument it can be shown that V (S(t)) ≤ c for all t ≥ 0. Hence, S(t) 2 is bounded for all t ≥ 0, and hence, the solution S(t) is bounded for all t ≥ 0.
Proof. The proof is similar to the proof of Theorem 5 using the function V : [0, ∞) × R n + → R given by and, hence, is omitted. (41) consisting of a netwok of six excitatory and three inhibitory neurons shown in Figure 4. The neural connectivity matrix A for this network is given by

Example 3 Consider the biological neural network given by Equation
Here, we assume that the time constant of all nine neurons are the same, and hence, L = 1 τ I 9 . Furthermore, we assume the functions f i (·), i = 1, 2, . . . , 9, are given by Equation (61) with the same maximum firing rate f max = 0.5.
Next, we construct A and L using Lemma 1 and solve the Linear Matrix Inequalities (LMIs) given by Equations (54) and (65) for the neuron time constants τ = 10, 1, and 0.1. We use the MATLAB toolbox YALMIP for solving the LMIs. For τ = 10, there is no feasible solution to Equations (54) and (65), and, as shown in Figure 5, the synaptic drive of the inhibitory neurons oscillate whereas the synaptic drive of the excitatory neurons converge to a steady-state. For τ = 1, Equations (54) and (65) are also not satisfied and as shown in Figure 6, the synaptic drive of all the neurons converge to different steady state values, which implies that the network is not synchronized. Finally, for τ = 0.1, Equations (54) and (65) are satisfied with P , Q, and R, given by  Hence, the biological neural network is asymptotically synchronized. See Figure 7.          Example 4 When patients lose consciousness some parts of the brain are still functional (e.g., cardiac function) whereas other parts are suppressed. This can be captured by biological neural network models that exhibit partial synchronization wherein part of the system's state in synchronized and the other parts fire at normal levels. In addition, 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 [82]. There is also a second biphasic surge in the EEG power as the patient emerges from unconsciousness. Models that predict the aforementioned characteristics are of great clinical importance in providing the phenomenological trends of the anesthetic cascade.
To demonstrate partial synchronization of the model presented in Section 2, consider the biological neural network given by Equation (41) consisting of a population of neurons with six excitatory neurons E 1 -E 6 and six inhibitory neurons I 1 -I 6 . The neural connectivity matrix A for this network is given by which implies that the three inhibitory neurons I 4 , I 5 , and I 6 do not receive any inhibitory inputs. Here we assume that the excitatory neurons have a time constant λ E = 0.01 s and the inhibitory neurons have a prolonged time constant λ I = 1 s. Furthermore, we assume v E th = v I th = 0.1V for all neurons, and the functions f i (·), i = 1, 2, . . . , 12, are given by Equation (61) with the same maximum firing rate f max = 0.5. Figure 8 shows that as λ I increases, the excitatory neurons that are coupled to inhibitory neurons all go to a zero synaptic drive, whereas the inhibitory neurons that themselves are not coupled to inhibitory neurons synchronize to some finite value.

Stochastic Multistability for a Mean Field Synaptic Drive Firing Neuronal Model
Since the neurocortex contains on the order of 10 11 neurons, each supporting up to 10 5 synaptic contacts, we extend our dynamical system framework to a stochastic setting. Specifically, a large population of spiking neurons can be reduced to a distribution function describing their probabilistic evolution; that is, a function that captures the distribution of neuronal states at a given time [83]. In this section, we develop a stochastic field theory as in [84] for capturing neural activity in order to analyze system multistability. In the next section, we extend the results of this section to additionally address time delay functional models [85] in order to account for time delay and memory effects in inhibitory and excitatory networks.
In Sections 4 and 5 we developed deterministic multistability theory to explain the underlying mechanism of action for anesthesia and consciousness using a synaptic drive firing model framework [3]. In this section, we extend these results further by demonstrating multistability in the mean when the coefficients of the neuronal connectivity matrix are random variables. Specifically, we use a stochastic multiplicative uncertainty model to include modeling of a priori uncertainty in the coefficients of the neuronal connectivity matrix by means of state-dependent noise. The philosophy of representing uncertain parameters by means of multiplicative white noise is motivated by the Maximum Entropy Principle of Jaynes [86,87] and statistical analysis [88].
Maximum entropy modeling is a form of stochastic modeling wherein stochastic integration is interpreted in the sense of Itô to provide a model for system parameter uncertainty. The use of stochastic theory to model system parameter uncertainty has been used within a modern information-theoretic interpretation of probability theory [86,87,89]. In particular, rather than regarding the probability of an event as an objective quantity such as the limiting frequency of outcomes of numerous repetitions, maximum entropy modeling adopts the view that the probability of an event is a subjective quantity which reflects the observer's certainty to a particular event occurring. This quantity corresponds to a measure of information. The validity of a stochastic model for a biological neural network does not rely on the existence of an ensemble model but rather in the interpretation that it expresses modeling certainty or uncertainty regarding the coefficients of the neuronal connectivity matrix. Hence, a stochastic multiplicative uncertainty model utilizes state-dependent Gaussian white noise to represent parameter uncertainty by defining a measure of ignorance, in terms of an information-theoretic entropy, and then determining the probability distribution which maximizes this measure subject to agreement with a given model.
To develop stochastic multistability synaptic drive model, consider for simplicity of exposition the simplified mean field synaptic drive model where the coefficients of Equations (24) and (25), with f i (·) given by Equation (10), are randomly disturbed. Specifically, we assume that the initial value S(0) [S E (0), S I (0)] T is deterministic and contained in the nonnegative orthant of the state space, and consider the stochastic differential mean field synaptic drive model given by where represents Brownian motion, that is, a Wiener process, ν ∈ R indicates the intensity of the Gaussian white noise dw(t), and Here, we assume that every entry of the matrices A and L of the mean dynamics given by Equations (24) and (25) For the statement of the results in this section, we require some additional notation and definitions. Specifically, let (Ω, F, P) be the probability space associated with Equation (70), where Ω denotes the sample space, F denotes a σ-algebra, and P defines a probability measure on the σ-algebra F, that is, P is a nonnegative countably additive set function on F such that P(Ω) = 1 [90]. Note that Equation (70) is a Markov process, and hence, there exists a filtration {F t } satisfying F τ ⊂ F t ⊂ F, 0 ≤ τ < t, such that {ω ∈ Ω : S(t) ∈ B} ∈ F t , t ≥ 0, for all Borel sets B ⊂ R 2 contained in the Borel σ-algebra B. Finally, spec(X) denotes the spectrum of the square matrix X including multiplicity and R(Y ) denotes the range space of the matrix Y .
When every component of the vector AS(t)(dt + νdw(t)), t ≥ 0, is nonnegative, the stochastic dynamical system given by Equation (70) can be written as whereÃ A − L andÃ s ν(A − L). The multiplicative white noise model given by Equation (71) can be regarded as a parameter uncertainty model where dw(t) corresponds to an uncertain parameter whose pattern and magnitude are given byÃ s / Ã and Ã s , respectively. Note that if rank(A − L) < 2, then every point α ∈ N (Ã) is an equilibrium point of Equation (71). With a slight abuse of notation, we use E to denote the equilibrium set of Equation (70) or Equation (71). First, motivated by the definition of stochastic Lyapunov stability in [90], we have the following definition of stochastic semistability. For the statement of the next result, define dist(x, E) inf y∈E x − y . For a similar definition of stochastic semistability, see [91]. (70) is stochastically semistable if the following statements hold.

Definition 10 An equilibrium solution S(t) ≡ α ∈ E of Equation
The dynamical system given by Equation (70) is stochastically semistable if every equilibrium solution of Equation (70) is stochastically semistable. Finally, the system given by Equation (70) is globally stochastically semistable if it is stochastically semistable and P lim t→∞ dist(S(t), E) = 0 = 1 for every initial condition S(0) ∈ R 2 . If, alternatively, S(t) ≡ α ∈ E only satisfies i), then the equilibrium solution S(t) ≡ α ∈ E of Equation (70) is stochastically Lyapunov stable.
Definition 10 is a stability notion for the stochastic dynamical system given by Equation (70) having a continuum of equilibria and is a generalization of the notion of semistability from deterministic dynamical systems [63,92] to stochastic dynamical systems. It is noted in [92] that existing methods for analyzing the stability of deterministic dynamical systems with isolated equilibria cannot be used for deterministic dynamical systems with nonisolated equilibria due to the connectedness property of equilibrium sets. Hence, Definition 10 is essential for analyzing the stability of stochastic dynamical systems with nonisolated equilibria. Note that (i) in Definition 10 implies stochastic Lyapunov stability of an equilibrium, whereas (ii) implies almost sure convergence of trajectories to the equilibrium manifold.
Next, we extend the notion of multistability for deterministic dynamical systems defined in Section 4 to that of stochastic multistability for stochastic dynamical systems.
Stochastic multistability is a global stability notion for the stochastic dynamical system given by Equation (70) having isolated equilibria and/or a continuum of equilibria, whereas stochastic semistability is a local stability notion for the stochastic dynamical system given by Equation (70) having a continuum of equilibria. Hence, stochastic multistability is a stronger notion than stochastic semistability. The next result states a relationship between stochastic multistability and global stochastic semistability.
Next, recall from [93] that a matrixÃ ∈ R n×n is semistable if and only if lim t→∞ eÃ t exists. In other words,Ã is semistable if and only if for every λ ∈ spec(Ã), λ = 0 or Re λ < 0 and if λ = 0, then 0 is semisimple. Furthermore, ifÃ is semistable, then the index ofÃ is zero or one, and hence,Ã is group invertible. The group inverseÃ # ofÃ is a special case of the Drazin inverseÃ D in the case whereÃ has index zero or one [93]. In this case, lim t→∞ eÃ t = I q −ÃÃ # [93].
Lemma 3 LetÃ ∈ R n×n . If there exist n × n matrices P = P T ≥ 0 and R = R T ≥ 0, and a nonnegative integer k such that Proof. The proof is similar to that of Lemma 4.5 of [96] and, hence, is omitted.

Proof. By Proposition 3 it suffices to show thatÃ
A − L is semistable. Consider the deterministic dynamical system given byẋ where x(t) ∈ R 2 . Note thatÃ is semistable if and only if Equation (82) is semistable [93], and hence, it suffices to show that Equation (82)  Next, consider the nonnegative function If V(x) = 0 for some x ∈ R 2 , then PÃ k x = 0 and Lx = 0. Now, it follows from Lemma 3 that x ∈ N (Ã), whereas Lx = 0 implies x ∈ R(Ã), and hence, V(x) = 0 only if x = 0. Hence, V(·) is positive definite. Next, since LÃ =Ã −ÃÃ #Ã = 0, it follows that the time derivative along the trajectories of Equation (82) is given bẏ Finally, let x e ∈ N (Ã) be an equilibrium point of Equation (82) and consider the Lyapunov function candidate U(x) = V(x − x e ), which is positive definite with respect to x e . Then it follows that the Lyapunov derivative along the trajectories of Equation (82) is given bẏ Thus, it follows that x e is Lyapunov stable. Now, it follows from Theorem 3.1 of [63] that Equation (82) is semistable, that is, A is semistable. Finally, it follows from Proposition 3 that Equation (71) is stochastically multistable.

Example 5
In this example, we illustrate the stochastic multistability properties of the two-state nonlinear synaptic drive neuronal firing model given by Equation (70). Specifically, consider the mean field synaptic drive model given by Equation (70) n I A II = 0 V, and λ E = 10 ms, and let λ I vary. In this case, the system matrices in Equation (71) are given by Figure 9 shows the eigenvalues of A − L as a function of λ I . Note that for λ I < 0.9 ms or λ I > 0.93 ms, (A − L) is unstable, whereas for 0.9 ms < λ I < 0.93 ms, (A − L) is asymptotically stable. Clearly, rank (A − L) < 2 for λ I = 0.9 ms. Hence, it follows from Theorem 7 that the stochastic dynamical system given by Equation (71)  For our simulation, we use the initial condition S(0) = [0.1, 0.5] T . Figures 10 and 11 show the time response for the average excitatory and inhibitory synaptic drives, and the phase portrait for λ I = 0.9 ms with ν = 0.2. Furthermore, for λ I = 0.9 ms with ν = 0.2, Figure 12 shows a histogram of the limit points of log e S E (or, equivalently, log e 0.9S I ) over 10, 000 samples. Note that the mean and variance of log e S E is −1.6135 and 0.3152, respectively. Similar plots shown in Figures 10 and 11 are shown for λ I = 0.78 ms and λ I = 1.20 ms with ν = 0.2 in Figures 13-16. Finally, Figures 17 and 18 show similar simulations for the case where λ I = 0.9 ms ( i.e., (A − L) is semistable) and ν = 1. However, note that in this case the condition given by Equation (77) of Proposition 3 is not satisfied, and hence, the model exhibits instability.

S I
The trajectories of Equations (24) and (25) can exhibit unstable and multistable behaviors for different values of the parameters, which are similar to the simulation results for Equation (70).
Moreover, its averaging dynamics will be analogous to the results of the deterministic model given in [98].

A Synaptic Drive Firing Model with Time-Varying Delays and Stochastic Multiplicative Uncertainty
In this section, we extend the synaptic drive model developed in Section 2 to investigate the conditions that would lead to synchronization or neutral oscillators. In particular, we extend the biological neural network model of Section 2 to include time-varying delays and stochastic input uncertainty. The system uncertainty model involves a Markov process wherein stochastic integration is interpreted in the sense of Itô.
For the statement of the results of this section, we require some additional notation and definitions. Specifically, C([−τ, 0], R n ) with τ > 0 denotes a Banach space of continuous vector-valued functions mapping the interval [−τ, 0] into R n with topology of uniform convergence and designated operator norm given by |||ψ||| = sup −τ ≤θ≤0 ψ(θ) for ψ ∈ C([−τ, 0], R n ). Furthermore, let S t ∈ C((−∞, +∞), R n ) defined by S t (θ) S(t + θ), θ ∈ (−∞, 0], t ≥ 0, denote an (infinite dimensional) state at time t corresponding to the piece of trajectories S between −∞ and t, and assume v thi (t) ≡ 0. To capture communication delays in our biological neural network model given by Equation (9), where δ ij (t) denotes the continuous, time-varying time delay of the transmission signal from the jth neuron to the ith neuron at time t, δ ij (t) ≥ 0, t ≥ 0, and S j (t) denotes the jth component of S(t). The system delays δ ij (t) correspond to the times of the spike hitting the synapse and t is the time after the spike, and hence, these delays account for the distance traveled by the voltage spikes down the axon. We modify the biological neural network system given by Equation (9) to include the effects of stochastic perturbations as well as time delays. Specifically, we consider the model where φ(·) ∈ C C((−∞, 0], R n ) is a continuous vector-valued function specifying the initial state of the system given by Equation (84), w(t) = [w 1 (t), w 2 (t), . . . , w n (t)] T captures noise in the input voltage and is represented by Brownian motion, that is, an n-dimensional mutually independent standard Wiener process, and σ(S) = diag[σ 1 (S), σ 2 (S), . . . , σ n (S)] represents the state-dependent noise intensity matrix for the Gaussian white noise process dw(t). Henceforth, we consider Equation (84) as the model of the perturbed biological neural network.
Next, sinceŜ(t) defined by Equation (83) contains n(n−1) terms with different time delays, each term can be written as the product of an n × n-dimensional matrix and an n-dimensional vector. Specifically, for i = 1, 2, . . . , n, j = 1, 2, . . . , n, i = j, define i i (n−1)+j, i > j, and i i (n−1)+j−1, i < j, where i = 1, 2, . . . , n(n − 1), define δ i (t) δ i j (t), and define the matrix A i ∈ R n×n whose (i , j)th entry is A i j and all the other entries are 0. Thus, the ith term in Equation (83) can be replaced by A i S(t − δ i (t)), i ∈ {1, 2, . . . , n(n − 1)}. Hence, setting N = n(n − 1),Ŝ(t) can be written aŝ For the statement of the results in this section, we define the infinitesimal operator L : [0, ∞) × C((−∞, 0], R n ) → R associated with the stochastic process given by Equation (84), acting on the functional V : R × C → R, by For a two-times continuously differentiable function V : [0, ∞) × R n → R of the random variable S, the infinitesimal operator LV (t, S) is defined as [90] LV (t, S) lim where V (t, S) denotes the Fréchet derivative of V and V (t, S) denotes the Hessian matrix of V with respect to S at (t, S). The following lemma provides an explicit formula for the infinitesimal operator on two kinds of functionals using the ideas from Lemma 3.1 of [99].
Lemma 4 Consider the biological neural network given by Equation (84) and let where Then, the infinitesimal operator acting on V 1 : [0, ∞) × C → R and V 2 : [0, ∞) × C → R are given by Proof. For sufficiently small h > 0, where O(h) denotes higher-order terms in h. Substituting Equation (92) into Equation (86) yields Equation (90). The proof of Equation (91) is similar to the proof of Equation (90) and, hence, is omitted.
To develop a global synchronization property for the biological neural network system given by Equation (84), we introduce the notion of stochastic synchronization.
Here, we focus on mean-square synchronization.

Synchronization of Stochastic Biological Neural Networks
In this section, we develop sufficient conditions for global mean-square synchronization for the biological neural network given by Equation (84) with differentiable time delays using Barbalat's lemma and linear matrix inequalities (LMIs). Here, we assume that the noise intensity matrix function σ(S) has a linear growth rate; that is, there exists r > 0 such that where M is defined in Equation (60).
The following theorem provides sufficient conditions for global mean-square asymptotic synchronization of the biological neural network system given by Equation (84).
Proof. Consider the functional V : [0, ∞) × C → R given by It follows from Equation (87) and Lemma 4 that the infinitesimal operator LV (t, S t ) associated with the stochastic process given by Equation (84) is given by where andŜ(t) ∈ R n is defined by Equation (85).
Proof. The proof is similar to the proof of Theorem 8 using the functional V : [0, ∞) × C → R given by and, hence, is omitted.
The following corollary to Theorem 9 is immediate.
Remark 2 Note that Theorem 8 does not require that the time delays be bounded, whereas Theorem 9 and Corollary 1 hold for the case where the time delays are bounded.

Remark 3
It is important to note that if f i (·), i = 1, 2, . . . , n, in Equation (84) is replaced by Equation (10), then the results of Theorems 8 and 9 as well as Corollary 1 still hold.

Thermodynamics, Neuroscience, Consciousness, and the Entropic Arrow of Time
In this and the next section, we present some qualitative insights from the fields of thermodynamics and electromagnetic field theory that can potentially be useful in developing mechanistic models [100] that can explain the underlying mechanisms of action for general anesthesia and consciousness. Specifically, by merging the two universalisms of thermodynamics and dynamical systems theory with neuroscience, one can provide key insights into the theoretical foundation for understanding the network properties of the brain by rigorously addressing large-scale interconnected biological neuronal network models that govern the neuroelectric behavior of biological excitatory and inhibitory neuronal networks. In addition, electrical signals in the brain can generate electromagnetic fields that can cause a shielding effect between the thalamus and frontal cortex, which in turn can lead to the emergence of unconsciousness. Both these paradigms are the subject of ongoing research.
In current clinical practice of general anesthesia, potent drugs are administered which profoundly influence levels of consciousness and vital respiratory (ventilation and oxygenation) and cardiovascular (heart rate, blood pressure, and cardiac output) functions. These variation patterns of the physiologic parameters (i.e., ventilation, oxygenation, heart rate variability, blood pressure, and cardiac output) and their alteration with levels of consciousness can provide scale-invariant fractal temporal structures to characterize the degree of consciousness [101] in sedated patients. Here, we hypothesize that the degree of consciousness reflects the adaptability of the central nervous system and is proportional to the maximum work output under a fully conscious state divided by the work output of a given anesthetized state. A reduction in maximum work output (and cerebral oxygen consumption) or elevation in the anesthetized work output (or cerebral oxygen consumption) will thus reduce the degree of consciousness. Hence, the fractal nature (i.e., complexity) of conscious variability is a self-organizing emergent property of the large-scale interconnected biological neuronal network since it enables the central nervous system to maximize entropy production and dissipate energy gradients. Within the context of aging and complexity in acute illnesses, variation of physiologic parameters and their relationship to system complexity and system thermodynamics have been explored in [102][103][104][105][106][107].
Complex dynamical systems involving self-organizing components forming spatio-temporal evolving structures that exhibit a hierarchy of emergent system properties form the underpinning of the central nervous system. These complex dynamical systems are ubiquitous in nature and are not limited to the central nervous system. Such systems include, for example, biological systems, immune systems, ecological systems, quantum particle systems, chemical reaction systems, economic systems, cellular systems, and galaxies, to cite but a few examples. The connection between the local subsystem interactions and the globally complex system behavior is often elusive. These systems are known as dissipative systems [108,109] and consume energy and matter while maintaining their stable structure by dissipating entropy to the environment.
In the central nervous system billions of neurons interact to form self-organizing dissipative nonequilibrium structures [108][109][110]. The fundamental common phenomenon among these systems are that they evolve in accordance to the laws of (nonequilibrium) thermodynamics, which are among the most firmly established laws of nature. System thermodynamics, in the sense of [110], involves open interconnected dynamical systems that exchange matter and energy with their environment in accordance with the first law (conservation of energy) and the second law (nonconservation of entropy) of thermodynamics. Self-organization can spontaneously occur in such systems by invoking the two fundamental axioms of the science of heat. Namely, (i) if the energies in the connected subsystems of an interconnected system are equal, then energy exchange between these subsystems is not possible, and (ii) energy flows from more energetic subsystems to less energetic subsystems. These axioms establish the existence of a system entropy function as well as equipartition of energy [110][111][112] in system thermodynamics and synchronization [73] in neuronal networks; an emergent behavior in thermodynamic systems as well as neuroscience. Hence, in complex interconnected dynamical systems, self-organization is not a property of the systems parts but rather emerges as a result of the nonlinear subsystem interactions.
In recent research the authors in [110][111][112][113][114] combined the two universalisms of thermodynamics and dynamical systems theory under a single umbrella to develop a dynamical system formalism for classical thermodynamics so as to harmonize it with classical mechanics. While it seems impossible to reduce thermodynamics to a mechanistic world picture due to microscopic reversibility and Poincaré recurrence [110,115], the system thermodynamic formulation in [110,112] provides a harmonization of classical thermodynamics with classical mechanics. In particular, our dynamical system formalism captures all of the key aspects of thermodynamics, including its fundamental laws, while providing a mathematically rigorous formulation for thermodynamical systems out of equilibrium by unifying the theory of heat transfer with that of classical thermodynamics. In addition, the concept of entropy for a nonequilibrium state of a dynamical process is defined, and its global existence and uniqueness is established. This state space formalism of thermodynamics shows that the behavior of heat, as described by the conservation equations of thermal transport and as described by classical thermodynamics, can be derived from the same basic principles and is part of the same scientific discipline. Connections between irreversibility, the second law of thermodynamics, and the entropic arrow of time are also established in [110,112,113].
Building on the results of this paper, we propose to merge system thermodynamics with neuroscience to develop a theoretical foundation for the mechanisms of action of general anesthesia using the network properties of the brain by rigorously addressing the large-scale interconnected biological neuronal network model given by Equations (7) and (8). Even though simplified mean field models of the form given by Equations (24) and (25) have been extensively used in mathematical neuroscience literature to describe large neural populations, the complex large-scale interconnected system given by Equations (7) and (8) is essential in indentifying the mechanisms of action for general anesthesia. We postulate that unconsciousness is associated with reduced physiologic parameter variability, reflecting the inability of the central nervous system to adapt. The degree of consciousness is a function of the numerous coupling in the network properties of the brain that form a complex large-scale, interconnected system. Complexity here refers to the quality of a system wherein interacting subsystems self-organize to form hierarchical evolving structures exhibiting emergent system properties, and hence, a complex dynamical system is a system that is greater than the sum of its subsystems or parts. This complex system-involving numerous nonlinear dynamical subsystem interactions making up the system-has inherent emergent properties that depend on the integrity of the entire dynamical system and not merely a mean field simplified reduced-order model.
As in thermodynamics, neuroscience is a theory of large-scale systems wherein graph theory [116] can be used in capturing the (possibly dynamic) connectivity properties of network interconnections, with neurons represented by nodes, synapses represented by edges or arcs, and synaptic efficacy captured by edge weighting giving rise to a weighted adjacency matrix governing the underlying directed dynamic graph network topology. However, unlike thermodynamics, wherein energy spontaneously flows from a state of higher temperature to a state of lower temperature, neuron membrane potential variations occur due to ion species exchanges which evolve from regions of higher chemical potentials to regions of lower chemical potentials (i.e., Gibbs' chemical potential [114]). And this evolution does not occur spontaneously but rather requires a hierarchical (i.e., hybrid) continuous-discrete architecture for the opening and closing of specific gates within specific ion channels.
Merging our proposed dynamical neuroscience framework developed in Sections 2 and 3 with system thermodynamics [110,112,114] by embedding thermodynamic state notions (i.e., entropy, energy, free energy, chemical potential, etc.) within our dynamical system framework can allow us to directly address the otherwise mathematically complex and computationally prohibitive large-scale dynamical model given by Equations (7) and (8). In particular, a thermodynamically consistent neuroscience model would emulate the clinically observed self-organizing, spatio-temporal fractal structures that optimally dissipate energy and optimize entropy production in thalamocortical circuits of fully conscious patients. This thermodynamically consistent neuroscience framework can provide the necessary tools involving multistability, synaptic drive equipartitioning (i.e., synchronization across time scales), energy dispersal, and entropy production for connecting biophysical findings to psychophysical phenomena for general anesthesia. In particular, we hypothsize that as the model dynamics transition to an anesthetized state, the system will involve a reduction in system complexity-defined as a reduction in the degree of irregularity across time scales-exhibiting partial synchronization of neural oscillators (i.e., thermodynamic energy equipartitioning). This would result in a decrease in system energy consumption (myocardial depression, respiratory depression, hypoxia, ischemia, hypothension, venodialtion), and hence, a decrease in the rate of entropy production. In other words, unconsciousness is characterized by system decomplexification, which is manifested in the failure to develop efficient mechanisms to dissipate energy thereby pathologically retaining higher internal (or local) entropy levels.
The human brain is isothermal and isobaric, that is, the temperatures of the subnetworks of the brain are equal and remain constant, and the pressure in each subnetwork also remains constant. The human brain network is also constantly supplied with a source of (Gibbs) free energy provided by chemical nourishment of the blood to ensure adequate cerebral blood flow and oxygenation, which involves a blood concentration of oxygen to ensure proper brain function. Information-gathering channels of the blood also serve as a constant source of free energy for the brain. If these sources of free energy are degraded, then internal (local) entropy is produced.
In the transition to an anesthetic state, complex physiologic work cycles (cardiac respiratory pressure-volume loops, mithochondrial ATP production) necessary for homeostasis follow regressive diffusion and chemical reaction paths that degrade energy production and decrease the rate of entropy production. Hence, in an isolated large-scale network (i.e., a network with no energy exchange between the internal and external environment) all the energy, though always conserved, will eventually be degraded to the point where it cannot produce any useful work (oxygenation, ventilation, heart rate stability, organ function). Hence, all motion would cease leading the brain network to a state of unconsciousness (semistability) wherein all or partial [117] subnetworks will possess identical energies (energy equipartition or synchronization) and, hence, internal entropy will approach a local maximum or, more precisely, a saddle surface in the state space of the process state variables. Thus, the transition to a state of anesthetic unconsciousness involves an evolution from an initial state of high (external) entropy production (consciousness) to a temporary saddle state characterized by a series of fluctuations corresponding to a state of significantly reduced (external) entropy production (unconsciousness).
In contrast, in a healthy conscious human entropy dissipation occurs spontaneously and, in accordance with Jaynes' maximum entropy production principle [86,118], energy dispersal is optimal leading to a maximum entropy production rate. Hence, low entropy in (healthy) human brain networks is synonymous with consciousness and the creation of order (negative entropy) reflects a rich fractal spatio-temporal variability which, since the brain controls key physiological processes such as ventilation and cardiovascular function, is critical for delivering oxygen and anabolic substances as well as clearing the products of catabolism to maintain healthy organ function. And in accordance with the second law of thermodynamics, the creation and maintenance of consciousness (internal order-negative entropy) is balanced by the external production of a greater degree of (positive) entropy. This is consistent with the maxim of the second law of thermodynamics as well as the writings of Kelvin [119], Gibbs [120,121], and Schrödinger [122] in that the creation of a certain degree of life and order in the universe is inevitably coupled by an even greater degree of death and disorder [110].
In a network thermodynamic model of the human brain, consciousness can be equated to the brain dynamic states corresponding to a low internal system entropy. In recent research [112], the author shows that the second law of thermodynamics provides a physical foundation for the arrow of time.
In particular, the author shows that the existence of a global strictly increasing entropy function on every nontrivial network thermodynamic system trajectory establishes the existence of a completely ordered time set that has a topological structure involving a closed set homeomorphic to real line, which establishes the emergence of the direction of time flow. Thus, awareness of the passage of time is a direct consequence of regressive changes in the continuous rate of entropy production taking place in the brain and eventually leading (in finite time) to a state of no entropy production (i.e., death). Since these universal regressive changes in the rate of entropy production are spontaneous, continuous, and decreasing in time, human experience perceives time flow as unidirectional. However, since the rate of time flow and the rate of system entropy regression (i.e., free energy consumption) are bijective (i.e., one-to-one and onto), the human experience of time flow is subjective.
During the transition to an anesthetized state, the external and internal free energy sources are substantially reduced or completely severed in part of the brain leading the human brain network to a semistable state corresponding to a state of local saddle (stationary) entropy. Since all motion in the state variables (synaptic drives) ceases in this unconscious (synchronized) state, our index for the passage of time vanishes until the anesthetic wears off allowing for an immediate increase of the flow of free energy back into the brain and other parts of the body. This, in turn, gives rise to a state of consciousness wherein system entropy production is spontaneously resumed. Merging system thermodynamics with multistability theory and mathematical neuroscience with the goal of providing a mathematical framework for describing the anesthetic cascade mechanism is the subject of current research. In addition, connections between thermodynamics, neuroscience, and the arrow of time [110,112,113] are also being explored to develop an understanding on how the arrow of time is built into the very fabric of our conscious brain.

An Electromagnetic Field Theory of Consciousness
The cerebral cortex can be modeled by a columnar topology, where thousands of 0.3-0.6 mm-wide cortical macrocolumns are compactly organized in a parallel configuration to create the approximately 4 mm-thick cerebral cortex [123]. Macrocolumns are bundles of approximately 100, 000 neurons which are arranged in six distinguishable layers and create a six-layered structure of the cortex; see Figure 21. Layers I, II, III, V, and VI are mainly composed of excitatory pyramidal cells, whereas Layer IV is almost free of pyramidal cells and is mainly composed of stellate cells. The apical dendrites of pyramidal cells perpendicularly extend towards the cortex superficial layer (Layer I), where they receive inputs from other parts of the brain. Schematic representing the connective topology within the six-layer cortical macrocolumn. Information is communicated inside the brain by electrical signals, or more specifically, by the creation and transmission of ions. The electrical signal received by dendrites at synapses travels along the dendrites to the cell body, and the resulting cell-generated signal travels to another cell along the axon. Thus, the spatial topology of the cortex dictates the paths along which the ions can flow through and the locations where they can reside. Based on the specific columnar topology of the cerebral cortex described above, the major portion of the electric current inside the cortex, which is conducted by pyramidal cells, flows radially through the cortex.
The almost evenly distributed short-range current flow through stellate cells in Layer IV do not significantly contribute to the global orientation of the current flow through the cortex. When apical dendrites of pyramidal cells receive excitatory inputs through synapses in Layer I, positive ions flow through them from synaptic clefts toward cell bodies, creating a deficit of positive ions at the synaptic clefts. The positive ions flowing through apical dendrites generate a transient radial current and, due to the capacitive characteristics of the cell membrane, accumulate around the cell membrane. As a result, macroscopic electrical activities in cerebral cortex possess four global behaviors (see Figures 21 and 22). Namely, (i) at Layer I, the potential transition due to the activity of neurons is negative; (ii) at layers II, III, V, and VI, the potential transition due to the activity of neurons is positive; (iii) no significant potential transition is observed in Layer IV; and (iv) transient currents due to neuron activity flows radially inside the cortex. Sparse activities of neurons in a region of the cortex do not generate significant electromagnetic effect in the cortex. However, when a large portion of neurons in a column or adjacent columns are hyperactivated, that is, fire synchronously, the resulting electromagnetic field due to the structured movement and ion configuration described above is strong enough to create significant electromagnetic interference inside the cortex; see Figure 22.
The state of unconsciousness due to the induction of general anesthesia can potentially be due to the strong electromagnetic pattern formation in the cerebral cortex. Available physiological evidence strongly suggests that the excitatory input signals received from the reticular activating system [123] are crucial in maintaining the conscious state of cerebral activities. Anesthetics are believed to prolong the post synaptic potential of inhibitory neurons, which mainly reside in Layer IV of the cortex [123]. Inhibitory neurons provide negative feedback signals for the population of neurons of cortical columns, and hence, can play an essential role in the stability of the highly interconnected network of neurons as well as the global behavior of these neurons. In particular, anesthetics can reduce the inflow of free energy to the brain leading to global or partial [117] synchronization of neuron firing.
Synchronization in the population of neurons present in Layers II, III, V, and VI of the cortex can create a strong electromagnetic pattern. The resulting electromagnetic field can affect the reception of signals from the reticular activating system by reducing the conduction velocity, partially or totally blocking, or even reflecting the incoming signals. Consequently, the cortex loses a large amount of its crucial background wash of input signals, which can result in loss of consciousness [15].
Even though electromagnetic field models have been proposed to explain consciousness (see [124] and the references therein), these models are largely based on emperical observation and lack precise mathematical formulations. A mathematical framework for fostering precision, completeness, and self-consistency in understanding the anesthetic cascade in the human brain using system thermodynamics and electromagnetic field theories are the subject of current research by the authors.

Conclusions
With advances in biochemistry, molecular biology, and neurochemistry there has been impressive progress in the understanding of the function of single neurons. Using the example of the mechanism of action of general anesthesia, the past decade has seen a remarkable explosion of our understanding of how anesthetic agents affect the properties of neurons. However, despite this advance, we still do not understand how molecular mechanisms translate into the induction of general anesthesia at the macroscopic level. In particular, 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. This clinical property leads to consideration of anesthesia as a nearly binary (on-or-off) variable, and the relationship between the concentration of an anesthetic agent in the central nervous system and the anesthetic state is described in terms of the probability of responsiveness as a function of anesthetic concentration [10]. In clinical studies, the typical observation is that at low concentrations of anesthetic agent the probability of responsiveness (to noxious stimuli) is high, possibly unity. Then as the anesthetic concentration increases there is a sharp transition to a probability of responsiveness that is low and possibly zero.
In this paper, we developed a synaptic drive firing rate model to model the central nervous system as a (possibly discontinuous) autonomous dynamical system and showed that the transition to the anesthetic state exhibits multistability; that is, the system exhibits multiple attracting equilibria under asymptotically slowly changing system parameters directly affected by anesthetic concentrations. In future research we plan to merge system thermodynamics, multistability theory, and dynamic network graph-theoretic notions to develop a framework for understanding central nervous system behavior characterized by abrupt transitions between mutually exclusive conscious and unconscious states.
Such phenomena are not limited to general anesthesia and can be seen in biochemical systems, ecosystems, gene regulation and cell replication, as well as numerous medical conditions (e.g., seizures, epilepsy, schizophrenia, hallucinations, etc.) and are obviously of great clinical importance but have been lacking rigorous theoretical frameworks. The primary impact of such frameworks will be to allow for the development of models that go beyond words to dynamic equations, leading to mathematical models with greater precision and self-consistency. Mathematical formulations enforce self-consistency and while "self-consistency is not necessarily truth, self-inconsistency is certainly falsehood." And within the context of general anesthesia, a dynamical system formulation for neuroscience can foster the development of new frameworks that will allow us to interpret experimental and clinical results, connect biophysical findings to psychophysical phenomena, explore new hypothesis based on the cognitive neuroscience of consciousness and develop new assertions, and improve the reliability of general anesthesia. In addition, such a framework can establish a scientific basis for new metrics of anesthetic depth by making the assessment of consciousness a mechanistically grounded tool.