Hysteresis in anesthesia and recovery: Experimental observation and dynamical mechanism

The dynamical mechanism underlying the processes of anesthesia-induced loss of consciousness and recovery is key to gaining insights into the working of the nervous system. Previous experiments revealed an asymmetry between neural signals during the anesthesia and recovery processes. Here we obtain experimental evidence for the hysteresis loop and articulate the dynamical mechanism based on percolation on multilayer complex networks with self-similarity. Model analysis reveals that, during anesthesia, the network is able to maintain its neural pathways despite the loss of a substantial fraction of the edges. A predictive and potentially testable result is that, in the forward process of anesthesia, the average shortest path and the clustering coefficient of the neural network are markedly smaller than those associated with the recovery process. This suggests that the network strives to maintain certain neurological functions by adapting to a relatively more compact structure in response to anesthesia.

The dynamical mechanism underlying the processes of anesthesia-induced loss of consciousness and recovery is key to gaining insights into the working of the nervous system. Previous experiments revealed an asymmetry between neural signals during the anesthesia and recovery processes. Here we obtain experimental evidence for the hysteresis loop and articulate the dynamical mechanism based on percolation on multilayer complex networks with self-similarity. Model analysis reveals that, during anesthesia, the network is able to maintain its neural pathways despite the loss of a substantial fraction of the edges. A predictive and potentially testable result is that, in the forward process of anesthesia, the average shortest path and the clustering coefficient of the neural network are markedly smaller than those associated with the recovery process. This suggests that the network strives to maintain certain neurological functions by adapting to a relatively more compact structure in response to anesthesia.

I. INTRODUCTION
The origin of consciousness is an unsolved mystery in nature [1][2][3][4][5][6]. To understand the neural physics of consciousness remains a challenging problem, requiring interdisciplinary efforts among researchers from disciplines such as neuroscience, physics, nonlinear dynamics and complex systems. Anesthesia-induced loss of consciousness and the possible recovery open a door to probing into the neural, physical, and dynamical mechanisms of consciousness [7][8][9][10]. Previous experimental studies provided evidence for the existence of a hysteresis phenomenon underlying the dynamics of loss and recovery of consciousness induced by anesthesia [11][12][13][14]. For example, signatures of a hysteresis were observed in behavioral experiments on mice and drosophila subject to injection of anesthetic and it was found that pharmacokinetics alone were not sufficient to explain the emergence of the hysteresis [12]. The concept of neural inertia was then proposed [12]. Neural inertia and hysteresis in humans were investigated through behavior indices [15,16] such as moments of loss and recovery of responsiveness, as well as through EEG measurements [17], providing guidance to clinical anesthesia practices [18]. In the work suggest- * juewang@xjtu.edu.cn † huangzg@xjtu.edu.cn ing neural inertia and hysteresis in humans based on the slow-wave activity in EEG through the dose-response relationship associated with induction and emergence [17], saturation in the slow-wave activity during anesthesia was found to depend on the age of the subject. This indicated that neural inertia and hysteresis might have a neurobiological basis in terms of the number of neurons and the synaptic density which typically deteriorate with age.
In the anesthesiology literature, hysteresis is an established phenomenon describing the changes in the patient EEG patterns as a function of drug concentration, which can be explained by the pharmacokineticpharmacodynamic (PKPD) model to some level. A source of debate is the extent to which the hysteresis loop can (or cannot) be nulled out via appropriate adjustments to the PKPD model used to compensate for delays between the measurement site (either blood concentration for intravenous drugs such as propofol, or endtidal lung concentration for inhalational drugs such as isoflurane) and the effect site, i.e., the brain. There has been extensive work on direct observation of the hysteresis effect by Kuizenga et al. [19][20][21], which suggests that the PKPD model is likely an oversimplification of the complex processes involved in generating the hysteresis effect in human anesthesia. The hysteresis effect is thus likely to have additional causes. In parallel, there were modeling efforts to understand the experimental results.
For example, a theory based on phase transition was proposed by Steyn-Ross et al. [22][23][24], providing an initial explanation of both biphasic excitation and hysteresis from the perspective of neuronal populations. Later, a percolation model based on reduced stochastic dynamics was proposed to describe the mechanism of general anesthesia [25]. Quite recently, a multistate Markov chain model was introduced to interpret neural inertia [26].
The experimental evidence supporting the concept of neural inertia was based on relatively indirect measurements, e.g., behavior indices or EEG, of the neural system through noninvasive detection. From the perspective of modeling, there is still a lack of understanding of neural inertia and the corresponding hysteresis from a network perspective. The purpose of this paper is to address the two issues. In particular, we have carried out animal experiments and obtained evidence that a hysteresis loop arises unequivocally in the processes of anesthetic administration and recovery after excluding pharmacokinetic interference effects. Our result is based on the local field potential (LFP) data that represent direct measurement of the neural system. Our measurements have revealed that the depth of anesthesia is state dependent. Motivated by the fact that state dependence is common in complex dynamical systems [27][28][29] and by the existent theoretical framework of modeling general anesthesia as a first-order phase transition in the cortex [24] in which the hysteresis is associated with state dependence, we develop a complex-network based dynamical mechanism to probe into the origin of the hysteresis phenomenon. The class of networks we construct belongs to multilayer complex networks with self-similarity, which involve the interactions of hierarchical units in the neural system responsible for anesthesia and recovery. Computations reveal that, during anesthesia, the network is able to maintain its neural pathways despite the loss of a substantial fraction of the edges. A predictive and potentially experimentally testable result from our network model is that, for a given anesthesia level during the forward process, the corresponding characteristics of the network, such as the clustering coefficient and average degree, can be different from those in the recovery process. A biological implication is that the network strives to maintain certain neurological functions by adapting to a relatively more compact structure in response to dramatic external disturbances such as anesthesia -possibly a survival strategy naturally gained during evolution of the nervous system. While the hysteresis phenomenon arising in the anesthesia-recovery cycle can be explained at the population level of neurons, our paper presents an alternative approach to understanding the phenomenon from a network perspective.

II. EXPERIMENTAL PROCEDURE AND RESULTS
Our experiments were conducted on two male six to eight week old C57BL/6 mice at the time of surgery. Animals were housed in a standard environment on a 12/12-h light/dark cycle with light on at 07:00, and they were allowed ad libitum access to water and food. The use and care of animals followed the guidelines of the Xi'an Jiaotong University Animal Research Advisory Committee, and all efforts were made to minimize animal suffering.
Mice were anesthetized using isoflurane (4% for induction, 2% for maintenance in 100% O 2 during surgical procedures) and were placed in a stereotaxic apparatus. A small animal heating pad was used for the maintenance of body temperature at 37.0 ± 1.0 • C. Microelectrodes of impedance 0.5MΩ were implanted on the right medial prefrontal cortex (1.7mm anterior to bregma, 0.4mm lateral to the midline, and 1.5mm below the brain surface), which were used for LFP recordings. For all the recordings, two screws implanted in the occipital bone above the cerebellum were used as the ground and reference. The whole implant was fixed with glass ionomer dental cement.
After surgery, all animals were allowed to recover for three to six days before undergoing LFP recordings. LFPs were recorded extracellularly using a 128-channel data acquisition system (Cerebus, Blackrock Microsystems, Salt Lake City, USA). Animals were placed in the induction chamber and connected to the data recording cables, and oxygen was administered. LFP recordings started 10 min before the induction of anesthesia. Isoflurane concentrations were then stabilized at 1.5, 2.5, 1.5, and 0%, each for 20 min, to allow equilibration between inspired and end-tidal concentration and thus eliminate the pharmacokinetic interference. The time duration of 20 min is sufficient to ensure the equilibration between the inspired and end-tidal concentrations to eliminate the pharmacokinetics interference factors. The burst suppression ratio (BSR) values shown in F1(b) were obtained from a 5-min time interval after a 20-min transient period so as to avoid the effect of unstable drug density during the transient phase of pharmacokinetics. More specifically, the time series of LFPs used to calculate the BSR value in F1(b) were recorded when the effect of isoflurane to the neural system already approaches equilibrium after 20 min, as verified through the statistical behaviors of the signals. LFPs were collected at a sampling frequency of 1kHz, amplified (300×), and band-pass filtered (0.3-500Hz).
F1(a) shows the measured LFPs from one of the mice during anesthesia and recovery. In the rising and falling phases of anesthetic concentration, the BSR [30][31][32][33][34] of the LFPs is different at the same level of concentration, signifying a hysteresis, where BSR is defined as the fraction of the duration of the suppression state in the total time interval (a suppression state is defined as the LFP signal being between -0.1 and 0.1mV for a continuous time of at least 100ms [31]). As shown in F1(b), the LFP during the falling phase has a larger burst suppression ratio, suggesting that the mouse is still in deep anesthesia. The hysteresis phenomenon shown in F1 cannot be fully interpreted by pharmacokinetics, leading to the introduction of the concept neural inertia [12]. Measurements of the second mouse gave essentially the same result. Illustration of the proposed network model. The whole system has a hierarchical multilayer structure. Each node in a lower layer is connected to four nodes in the adjacent upper layer (as shown in the dashed box). Within each layer, the network can be scale free, small world, or random. The input signal is applied to the root node and the output signal is taken at a randomly selected node in the top layer.

A. Model construction
To interpret the experimental results, we develop a dynamical network model based on physical considerations of the neural network structure and the dynamics underlying the processes of anesthesia and recovery. For convenience, we call the process of applying anesthetic "forward" and that of recovery "backward." Structure-wise, in the two processes, many neural units are involved. Regarding each neural unit as a node and the connection strength between a pair of units as a weighted edge, we construct a class of multilayer brain neural networks [25]. An example is shown in F2, where the network consists of five layers with a self-similar, fractal-like structure. The first (bottom) layer has only one node, the root node, that serves as the modulator of the brain (thalamus) and provides the driving force for the whole brain. For the other layers, in terms of the possible structures of neuron connections, we consider general complex network structures such as the scale-free [35], small-world [36], or random [37] topology, to describe the nodal connections within each layer. The input stimulation signal is applied to the root node and the corresponding output signal is taken at a randomly selected node in the top layer.
We construct the network dynamics model based on the following intuitive reasoning. In a typical setting, in the forward process, as the amount of administrated anesthetic is increased, loss of consciousness occurs almost instantaneously at a certain moment rather than gradually. Likewise, during the backward (recovery) process, awakening tends to occur abruptly. The empirical observation suggests that both the forward and backward processes can be described by percolation dynamics on the underlying complex neural network [25,38,39], corresponding to progressive edge removal from and addition into the network, respectively. During the forward process, failures of connections between the neural units occur as a result of the action of anesthetic, depriving the units and those connected to them of the ability to transmit and integrate information. Intuitively, the connections among the more active nodes are more difficult to break than those among less active nodes, which correspond to the larger or smaller degrees of nodes from the aspects of topology. We thus assume that, at a certain anesthesia level, the connection failure between a pair of nodes is random with a probability determined by the degrees of both nodes. Effectively, there is mutual maintenance among the hub units to sustain the vital neural connections in the network. Neurological experiments indeed indicated that the normal release of neurons is associated with maintenance and consolidation of the synaptic connections and thus has a positive effect on the formation and stable survival of the entire neural network [40]. More active neurons with a higher firing rate typically have a large in-degree and thus are more likely to form and maintain their connections with other nodes [41]. Similarly, during the backward process, edges are gradually restored in the network. For a pair of nodes, the probability of recovery is determined by the current connectivity of the nodes. Consequently, the connections between large-degree nodes in the presently "broken" network are easier to recover but these nodes may not appear as the potential hub nodes, especially in the early stage of the recovery process.
The physical reasoning supporting our construction of the dynamical neural network suggests the following breaking and recovery probabilities associated with an edge of end nodes i and j: p b ij = 1 − (S i + S j )/(S max + S submax ) and p c ij = (S i + S j )/(S max + S submax ), respectively, where S i = k w ik characterizes the present strength of node i and S max and S submax are the maximum and the second largest strengths at the present moment in the network, respectively. The forward and reverse edge weights of the network are equal: w ij = w ji . The level of anesthetic can be characterized by the failure ratio ρ -the fraction of failed edges.

B. Emergence of hysteresis loops in various aspects of the dynamical neural network
Our main point is that the experimentally observed hysteresis in F1 can be understood as the result of simultaneous emergence of hysteresis loops in the characterizing quantities of the dynamical network, which is es-  FIG. 3. Emergence of a hysteresis loop in the network structure associated with the forward and backward processes. The hysteresis is revealed by the Kullback-Leibler divergence defined as the ratio of the value distribution of the output signal to that of the input signal vs the failure ratio ρ for (a) fixed edge weight w = 1 and (b) uniformly distributed random edge weight w ∼ U (0, 2). Each data point is the result of averaging over 100 statistically independent realizations. For a single anesthesia-recovery process, the transition from consciousness to unconsciousness is abrupt, and vice versa. The input is a Gaussian signal N(10 + sin (2πt/100), 1) with timevarying mean. Results from three types of network models are included: scale-free networks (SFNs), random networks (RNs), and small-world networks (SWNs). The abbreviations "Pro", "Net", "Ane" and "Rec" stand for process, network, anesthesia, and recovery, respectively. The same abbreviations are used in subsequent figures.
tablished through a systematic analysis of the dynamical network.
We first demonstrate the emergence of a hysteresis in the network structure. We use the informa-tion transmission capability to measure the awareness of the brain's nervous system, which can be characterized by the relative entropy, also known as Kullback-Leibler divergence [42], between the value distributions of the input and output signals: where p α (i) and p α (o) are the probabilities for the input and output signals to take on the value α, respectively, and α is calculated by the bin average. F3 shows the D KL (o||i) curves for the forward and backward processes. It can be seen that the two curves form a hysteresis loop enclosing a finite area, indicating that the mechanism of mutual maintenance between the network nodes can result in quite different network structures for the forward and backward processes at the same anesthesia level, especially at the intermediate levels. In general, different network structures lead to different information transmission capabilities and further to the experimentally observed hysteresis in the burst suppression ratio of LFPs as demonstrated in F1. The hysteresis in the network structure is robust, as it arises regardless of the topology of the network layers. Nonetheless, for randomly distributed edge weights, the area of the hysteresis loop is smaller than that for fixed weights.
We next demonstrate the emergence of another type of hysteresis, one manifested as the difference in the network connectivity between the forward and backward processes. F4(a) and F4(b) show E inter /E and E intra /E for the forward and backward processes versus the failure ratio ρ, where E inter , E intra , and E are the numbers of interlayer and intralayer edges, and all edges in the network, respectively. In the forward process, E intra decreases rapidly with the deepening of anesthesia since a large number of fragile edges with small degree ends are typically intralayer edges, while E inter decreases quite slowly. As a result, the ratio E inter /E increases monotonically with ρ while E intra /E decreases monotonically, as indicated by the dashed curves in F4(a) and F4(b), respectively. For the backward process where edges are gradually restored in the network, a pair of nodes with large degree values are more likely to be connected, and many such potential links are of the interlayer type to connect the hub nodes across the layers. As the value of ρ gradually decreases due to continuous reduction in the anesthesia level, E inter increases rapidly but E intra exhibits a slow increase at first. However, the potential interlayer edges are not many in comparison with the total number of possible edges in the network. When most of the interlayer edges have been restored, the continued increase in E inter tends to slow down, leading to a rapid increase in E intra . As a result, during the entire backward process, E inter /E first rises and then falls, as indicated by the solid curves in F4(a) (from right to left), while E intra /E first falls and then rises, as revealed by the solid curves in F4(b) (from right to left). Mutual maintenance between nodes causes the forward and backward processes to have different edge-connecting dynamics, leading to different network structures for the two processes at the same anesthesia level (medium ρ). During the entire backward process, the restoration of edges is comparatively random, and the small cliques formed by these restored edges are not well organized to support percolation between thalamus and cortex. The underlying dynamics in the backward process are thus not compatible with the requirement of consciousness recovering. The end result is a hysteresis that emerges in the generation of network pathways between the forward and backward processes. The hysteresis in the network structure manifests itself as hystereses in the characterizing quantities of the network, such as the average shortest path L from the input to the output node and the clustering coefficient C, as shown in F5(a) and F5(b), respectively. For intermediate values of the ρ, the values of L and C associated with the forward process are significantly lower than those with the backward process, indicating that, under the same anesthesia level, the underlying neural network of the forward process transmits information more rapidly with more straightforward and less redundant network connectivity than that with the backward process. Biologically, this may be interpreted as a kind of functional maintenance that the brain nervous system struggles to achieve in response to threats.

IV. DISCUSSION
To summarize, we have obtained experimental evidence of the emergence of a hysteresis loop underlying the anesthesia and recovery processes. To preclude possible interference from pharmacokinetics, we have carried out direct measurement of the local field potentials of the neural system through invasive detection. We have developed a physical understanding based on multilayer, self-similar networks with a complex topology. The specific neurophysiological mechanism of the hysteresis is unknown at the present, but our model suggests that the mutual maintenance between certain units of the nervous system may play an important role in generating the hysteresis. The mutual interactions, due to the heterogeneous connectivity of the nervous system, can result in a significant difference in the dynamical evolution of neural connections between the anesthesia and recovery processes. The complex brain neural network is the result of long term evolution and has embodied some kind of efficiency and adaptability with functional robustness in response to external disturbances or threats [43][44][45][46]. One of our findings is that, in the forward process of anesthesia, the average shortest path and the clustering coefficient of the neural network are markedly smaller than those associated with the recovery process (F5). This suggests that the network strives to maintain certain neurological functions by adapting to a relatively more compact structure in response to anesthesia -a kind of dramatic external disturbance. This represents a survival strategy naturally gained as the result of evolution of the nervous system. Additionally, our network model predicts that the structural characteristics of the underlying network, such as the clustering coefficient and average degree, can be quite different during the anesthesia and recovery processes. An experimental test of this prediction is called for.
There has been experimental work indicating that noise tends to attenuate the hysteresis associated with the anesthesia-recovery process [26], which can actually be explained with our network approach. In particular, hysteresis is a phenomenon in multistable dynamical systems going through state transitions. Noise can enhance the probability for a state transition to occur, thereby weakening the hysteresis effect. From the point of view of either neuronal populations or network, the anesthesiarecovery process can be viewed as the transitions of the nervous system between two stable states, so the underlying system is effectively bistable. In the forward process, the robust connections among the large degree nodes mean that the main pathway between the input and output nodes can always be maintained, regardless of perturbations. In this sense, noise has little effect on the forward process. However, noise can have a significant effect on the backward process, which can be seen by noting that, in the nervous system, noise typically manifests itself as some kind of remote excitatory stimulus that specifically enhances the excitability of the neurons and accordingly synaptic connections. In our network model, conceptually the role of noise in enhancing neuronal connections is equivalent to adding extra edges into the network. The extra edges significantly increase the probability of the emergence of pathways between the input and output nodes, advancing the occurrence of a fully connected network state. The global effect of noise is thus to weaken the hysteresis.