Characterization of Phase Transition in the Thalamocortical System during Anesthesia-Induced Loss of Consciousness

The thalamocortical system plays a key role in the breakdown or emergence of consciousness, providing bottom-up information delivery from sensory afferents and integrating top-down intracortical and thalamocortical reciprocal signaling. A fundamental and so far unanswered question for cognitive neuroscience remains whether the thalamocortical switch for consciousness works in a discontinuous manner or not. To unveil the nature of thalamocortical system phase transition in conjunction with consciousness transition, ketamine/xylazine was administered unobtrusively to ten mice under a forced working test with motion tracker, and field potentials in the sensory and motor-related cortex and thalamic nuclei were concomitantly collected. Sensory and motor-related thalamocortical networks were found to behave continuously at anesthesia induction and emergence, as evidenced by a sigmoidal response function with respect to anesthetic concentration. Hyperpolarizing and depolarizing susceptibility diverged, and a non-discrete change of transitional probability occurred at transitional regimes, which are hallmarks of continuous phase transition. The hyperpolarization curve as a function of anesthetic concentration demonstrated a hysteresis loop, with a significantly higher anesthetic level for transition to the down state compared to transition to the up state. Together, our findings concerning the nature of phase transition in the thalamocortical system during consciousness transition further elucidate the underlying basis for the ambiguous borderlines between conscious and unconscious brains. Moreover, our novel analysis method can be applied to systematic and quantitative handling of subjective concepts in cognitive neuroscience.


Introduction
Many studies have characterized electrophysiological differences between conscious and unconscious states such as vegetative state and coma [1], anesthesia-induced unconsciousness [2]; however, the dynamic processes of electrophysiological activity during consciousness transitions have been poorly characterized. As a prototypical example, anesthesia is used to correlate altered neuronal function or cortical communication with loss of consciousness [3,4]. Steyn-Ross et al. [5] formulated a theoretical model for the cerebral cortex response to general anesthesia and predicted a discontinuous phase transition of electrical fluctuations in the cortex at the point of transition into unconsciousness, which was characterized by a divergence in cortical synchrony and hysteresis in transitional paths. Their analogy of anesthetic phase transition to a thermodynamic phase change such as water freezing presumes that the anesthetic concentration is an external condition that drives the cortical state to an ordered phase, such as the unconscious state, by reducing neuronal excitability. Liley and Bojak [6] later adopted mean-field theory to correlate neuronal population behaviors with physiologically observable parameters (e.g., electroencephalogram (EEG)), and Molaee-Ardekani et al. [7] added voltage-dependent slow motion of ions to the above models, predicting the existence of a broad regime of fluctuation between two states, which is a property of continuous phase transition. These models describe the mechanisms underlying why and how anesthesia drives cortical neurons from the up state (depolarized mode) to the down state (hyperpolarized mode). However, there is no quantifiable variable representing the onset of consciousness transition suggested in freely behaving individual animal, as far as we know.
It has been shown that the resonance oscillations with large amplitude in the thalamocortical network are correlated with the functional states of cortical neurons [8,9]. A synchronous alteration of thalamocortical neuron up (depolarized) and down (hyperpolarized) states blocks afferent information and disassociates top-down integration. In the present study, the generation and breakdown of resonant oscillations in the thalamocortical network associated with anesthesia were investigated in sensory and motor thalamocortical circuit, with a novel analysis method quantifying an order parameter which is 0 and 1 in information transmission and blocking modes, respectively. Our primary purpose was to characterize thalamocortical network phase transition in the induction and emergence of anesthesia, with the aim of locating the critical points of transition to unconsciousness or consciousness in the thalamocortical network. Consciousness is a subjective term, and our strategy is to delimit the unconsciousness as a loss of ability to integrate perceptual and motor events together with memory to respond adequately to a given situation, which is referred to primary consciousness by Edelman [10]. The moments of behavioral change induced by anesthetic administration was pinpointed with motion tracker and three-compartment modeling was employed to estimate the ketamine concentration in the brain. The dynamic behavior of the sensorimotor thalamocortical network was characterized according to external conditions such as drug induction time and anesthetic concentration in the brain.

Surgical Procedures
The mouse model was used in research and all surgical and experimental procedures were conducted in accordance with the guidelines for the Institutional Animal Care and Use Committee, following Act 1992 of the Korea Lab Animal Care Regulations and associated guidelines.
Experiments were conducted on ten male C57BL/66129 F1 hybrid mice (8-10 weeks; body weight 19-25 g). Surgical procedures were performed as described in a previous report [11]. In brief, mice under ketamine-xylazine anesthesia (120 and 6 mg/kg, respectively, i.p.) were implanted with two tungsten electrodes (Teflon-insulated tungsten wire 76.2/114.3 mm in bare/ coated diameter, A-M Systems, Sequim, WA, USA) positioned at the ventral lateral (1.06 mm posterior, 1.1 mm lateral and 3.5 mm ventral to bregma) and ventral posteromedial nuclei (1.82 mm posterior, 1.5 mm lateral and 3.7 mm ventral to bregma) of the thalamus; two cortical screw electrodes (chrome-plated stainless steel 3 mm in length and 1 mm in diameter, Asia Bolt, Seoul, Korea) positioned at the primary motor (0.74 mm anterior, 1.5 mm lateral and 0 mm ventral to bregma) and primary somatosensory cortices (1.82 mm posterior, 3.0 mm lateral and 0 mm ventral to the bregma); and a ground electrode on the interparietal bone. All electrode coordinates followed the mouse brain atlas [12]. After experiments, histology was performed to verify the thalamic placement ( Fig. 1), and all the data with tip of electrode located out of the targeted thalamic nuclei were excluded in the analysis. For remote administration of the anesthetic agent during electrophysiological recording, polyethylene tubing (10 cm in length, 1520/86 mm in outer/inner diameter; PE 100, Intramedic polyethylene tubing, Clay Adams, Sparks, MD, USA) was inserted in the intraperitoneal cavity and sutured.

Electrophysiological Recordings under Forced Walking
Electrophysiological recording was performed under a forced walking scheme to prevent EEG artifacts induced by the experimenter's intervention during anesthetic drug injection and to obtain behavioral reference points of brain state transitions [11]. Prior to experiments, each mouse was habituated to a treadmill (LE8708, Panlab, Spain) in a suspended condition for approximately 15 min and in a walking condition for approximately 5 min. While the animal was walking on the treadmill, baseline EEG signals were recorded for 10 min, followed by unobtrusive administration of the ketamine-xylazine cocktail (120 and 6 mg/ kg, respectively) through the preinstalled injection tube. Recording continued until the animal awakened and recovered its movement. The motion of the animal was simultaneously monitored by an accelerometer-based motion sensor (MMA7260Q, Freescale Semiconductor Inc., Austin, TX, USA) attached to the animal's head and used to detect behavioral reference points of loss and recovery of motion referred to as t LOM and t ROM , respectively, with respect to the time of drug administration (t ADM ). The speed of the lane was maintained at 5 cm/s throughout the recording period, which induced minimal walking behavior for the animal equivalent to a spontaneously behaving condition. The sides of treadmill were surrounded with a wall to keep the animals on the treadmill. After loss of motion, the animals lied down above the running lane against the wall continuously receiving tactile stimulus generated by friction from the lane surface.
Simultaneously with motion signal acquisition, EEGs were acquired in a monopolar scheme referenced to the ground electrode by an analog amplifier (8-16C, Grass Technologies, West Warwick, RI, USA) and digitized with an analog-digital converter (Digidata 1440A, Molecular Devices, Sunnyvale, CA, USA) at a sampling frequency of 10 kHz. Before analysis, EEG signals were band-pass filtered with cut-off frequencies of 0.5 and 100 Hz (linear-phase Bessel filter), and then the sampling frequency was reduced to 200 Hz by moving average. Each signal was normalized by its average power in the frequency range of 90-100 Hz during the quiescent period.

Phase Locking Value as an Index of Phase Synchronization in Delta Band
Phase synchronization was calculated with phase locking value (PLV) defined as where Dw xy (t n ) is the phase difference between signal x and y at time t n and N is the number of data points in a 10-s analysis window. 1-4 Hz band-pass filtered EEG signals (zero-phase filtering with FIR filter of order 400) were Hilbert transformed to yield instantaneous phases and the PLVs between motor-related thalamocortical pair and somatosensory-related pair were computed.

Dimension Reduction with Aggregate Time Units
To reduce the complexity of the EEG signals, we introduced aggregate time epochs with variable intervals reflecting the period of one slow fluctuation. The aggregate time was defined as the point where more than three neuronal signals got into negative peaks simultaneously. The negative peaks of the signals were designated as the points satisfying the criterion that the Hilbert phase of the 1-4 Hz band-pass filtered signal (zero-phase filtering with FIR filter of order 400) equaled -p with 15% tolerance.

Concept of Order Parameter and Quantification of Information Blocking
When we assume each region represents a functional unit in terms of neuronal information transfer, the oscillation patterns in each epoch reflect functional states of the brain. To each state of a brain system functional unit, there corresponds a ''state vector'' in the description of its resonant oscillation. For example, suppose that a functional brain unit is perfectly resonant within the thalamocortical system, rarely responding to any perturbation and resultantly blocking signal transmission. We can designate the corresponding state vector as D1T when all the functional units in the brain system are synchronously in a resonant state, whereas D0T when none of the units experiences resonant oscillation.
The notation DyT i is used for the state of the thalamocortical system at the i th epoch to represent the brain state as for N different functional units, g j . An operator, A, acting on the state vector to evaluate brain state can be written as a superposition of independent operator A j acting on the state vector of each functional unit, i.e., For each operator, a state value for the j th functional unit at the i th epoch can be determined as A j Dg j T i~aj,i Dg j T i , where a j,i is 1 or 0 for ordered and disordered states, respectively. Here, an order parameter of interest, m i can be defined as with the value of m i ranging from 0 to 1. The signal operator was defined as where H is a Heaviside step function, h is a predetermined orderdisorder threshold, and R 2 j is a dissimilarity function. Dissimilarity was calculated by comparing the length-and amplitude-normalized oscillations in each epoch with a template of a representative epoch from deep anesthesia, which was drawn by normalizing and averaging up to 500 epochs sampled from 10 min after the drug administration. Dissimilarity function for the i th epoch, R 2 i , was calculated as where , and p and q are least-square fitting parameters minimizing R 2 i . In the present study, h was determined by k-means clustering of a sampled R 2 distribution from wakeful and anesthetized periods (k = 2). R 2 i value is conceptually similar to a z-score where larger absolute z-scores indicate larger deviation from the mean value. A pattern of the i th epoch is similar to the normalized template when R 2 i is close to zero, and it is more deviant from the template as R 2 i has larger value.
To estimate the readiness for transition between wakeful state and anesthesia, the susceptibility of the order parameter to the anesthetic concentration during anesthetic induction and recovery was calculated as x H~dm m H =dc and x D~dm m D =dc, respectively, wherem m H andm m D are sigmoidal fitting functions of m during anesthetic induction and recovery, respectively, and c is the anesthetic concentration estimated from mathematical modeling (see below section for mathematical modeling of anesthetic concentration). Anesthetic induction was defined as the period in which the estimated anesthetic concentration was in the increasing phase before reaching its maximum, and recovery as the period during which the anesthetic concentration decreased after the maximum. The shape of the sigmoidal fitting function followed m m(c)~x 1z exp ({yczz) , and the fitting coefficients x, y, and z were Phase Transition during Loss of Consciousness PLOS ONE | www.plosone.org estimated in a least squares manner. During the anesthetic induction period, the more easily hyperpolarized neuron groups were expected to have a higher level of susceptibility. Hence, susceptibility during induction was referred to as the ''hyperpolarizing susceptibility,'' x H . By contrast, susceptibility during recovery gauges the ability to resume information transfer by breaking the resonant oscillations; therefore, it was referred to as the ''depolarizing susceptibility,'' x D :

Probability of Transition between States
The probability of transition between states, d, in a 10-s observation window was defined as where n D0T?D1T and n D1T?D0T are the number of state transitions switching from a disordered (up state) to an ordered state (down state) and vice versa, and n epochs is the total number of epochs in the observation period, with an epoch defined by the time interval between adjunct aggregate time points. A zero value of d indicates that the thalamocortical system stays in a single state during the observation period, and positive unity indicates that the thalamocortical system continuously switches its state.

Mathematical Modeling of Anesthetic Concentration in the Brain
Prior to further investigating the characteristics of state transition, determining whether time is the appropriate variable to describe the brain state is a crucial question. Assuming susceptibility to anesthesia depends on the concentration of the anesthetic in the cerebral circulatory system [13], the concentration of ketamine in the cerebral circulatory system was considered the reference parameter reflecting the endogenous state of the brain. However, directly measuring this concentration is difficult in practice due to the limited access to the dense network of blood vessels in the brain in the electrode-implanted mice. Indirect measurement of ketamine concentration through the peripheral plasma also has limitations as it inevitably affects EEG recording, since the needle's sting may distort EEG signals or induce EEG artifacts during blood sampling. To circumvent these limitations, a mathematical model based on the behavioral reference points, i.e., t ADM , t LOM , and t ROM , which are believed to be directly proportional to the cerebral levels of ketamine, was used to predict the ketamine concentration in the brain.
The mathematical model consists of three compartments. The first (C 1 ) includes the initial dilution volume, which is composed of the rapidly perfused tissues and organs such as the brain. The second compartment (C 2 ) involves the tissues and organs of the body that are less well perfused. Finally, an additional compartment (C a ) that denotes the site of ketamine administration is incorporated into the model. The dynamics equations for the time evolution of the ketamine concentrations in the C a , C 1 , and C 2 compartments, designated by c a , c 1 , and c 2 , respectively, are determined by the difference between corresponding influx and efflux rates. These are expressed in terms of various specific rates, which describe elimination of ketamine (k a ) and transitions of ketamine (k 12 and k 21 ). The details of the model system are described in Information S1 with the model parameters in Table S1.
The resulting equation for the ketamine concentration in C 1 after administration is given by: where k a is the rate constant for drug absorption from the administration compartment, k 21 is the rate constant for drug transition from C 1 to C 2 , a is the rate constant for ketamine distribution process, b is the rate constant for ketamine elimination, F is the bioavailability, D is the loading dose, V 1 is the volume of C 1 , and t 0 is thedrug administration time. k a and a were determined to be inverse-proportional to the intervals between t ADM and t LOM , and between t LOM and t ROM , respectively, which were considered as proxies for anesthetic effects. The estimated ketamine concentration c 1 was normalized by its maximum value to yield values between 0 and 1. Later the normalized anesthetic concentration was denoted as c without a subscript.

Delayed Generation and Advanced Breakdown of Thalamocortical Slow Oscillation with Respect to Loss and Recovery of Motion
In Fig. 2A, 15-s sample signals of EEG and motion show qualitative differences of brain states in the vicinity of t ADM , t LOM , deep anesthesia, t ROM , and resumption of walking. The EEG signals in the middle of anesthesia period have regularly-repeating large amplitude oscillations in delta band range (1-4 Hz). In delta band, phase synchronization values for the motor-related thalamus-cortex pair and the somatosensory-related thalamus-cortex pair were increased and kept high during the anesthetic period (Fig. 2B).
The change of electrophysiological brain state and its saturation did not synchronously match the change in behavior response under anesthesia: the change of the brain state lagged behind t LOM during induction of anesthesia, and this temporal relationship was reversed during emergence, i.e., by advancing t ROM . Figure 2C shows the subject-averaged order parameter m over rescaled time, which is normalized against the blackout time (i.e., period between t LOM and t ROM ). The induction time (i.e., period between t ADM and t LOM ) was 96.7636.8 s (mean6SD) and the blackout time was 33.268.3 min (mean6SD). As indicated by the response curve (Fig. 2C), a time delay occurred before the order parameter reached its saturation value. Although the saturation points were roughly estimated due to the fluctuation in their vicinity, our estimation of the time delay of the saturation point with respect to t LOM was 2.060.6 min, and the time advance of the saturation point with respect to t ROM was 8.664.8 min. The average duration of the steady state of the order parameter was 22.765.3 min.

Quick Induction and Slow Emergence of Anesthesia in the Time Domain
The brain state transition to the down state during induction of anesthesia was noted to occur significantly faster than the transition to the up state during emergence from anesthesia. The difference in the speed of transition is reflected by the absolute values of the time derivative m in Fig. 2D, which describes the speed of phase transition. The generation of global synchronization to resonant oscillations was found to occur responsively within  Figure 2E shows the transitional probability between up and down states, d, within a 10-s interval. Notably, the time of maximal transitional probability coincided with the time of maximal transitional speed during the induction period in Fig. 2D. By contrast, the transitional probability began to increase at the end of the steady state of m until t ROM and did not drop to zero even after t ROM . The large temporal fluctuation between up and down states implies that the direction of transition during the transitional period is not fixed, strongly suggesting the non-deterministic nature of consciousness transition around transitional period.

Path-dependent Reaction to Anesthesia Induction and Washout
As previously stated, we observed large individual variation in drug induction time and blackout time across ten animals. We assumed the drug induction and blackout times as the time constants for drug absorption and early wash-out, respectively, and estimated the ketamine concentration in the brain using a threecompartment model (Fig. 3A). The subject-averaged response functions of m as a function of anesthetic concentration c are plotted in Fig. 3B with curve fitting. The lag of m is observed in both transitional paths, showing a hysteresis loop. This observed path-dependent nature of brain systems during transitions in consciousness is consistent with theoretical predictions [5] and path-dependent behavioral responses to anesthetic dose [14]. Furthermore, it is interesting to note that the hyperpolarizing susceptibility, x H (i.e., Dm/Dc in the induction period) reached its maximum at the anesthetic concentration value of 0.87560.106, which is significantly higher than the critical value of 0.25760.044 for the depolarizing susceptibility, x D (i.e., Dm/Dc in the emergence period) as shown in Fig. 3C. This result implies that a certain level of anesthetic concentration is needed for global synchronization, but neuronal recruitment grows extensively near the point of highest anesthetic concentration. The lower threshold for depolarization compared to the threshold for hyperpolarization suggests that once neurons are hyperpolarized, they have a tendency to maintain their resonant oscillations.

Temporary Metastable State of Global Activity during Transitional Periods
The calculated probability distribution of m as anesthetic concentration changes is shown in Fig. 3D. Interestingly, mixed states exist near transitional points with fractional values of m, e.g., near 0.8 in the increasing phase of c and near 0.2 in the decreasing phase of c. These bimodal peaks observed at the transitional periods imply the existence of spatially and temporally local up or down states in the thalamocortical network.
A close look at the microscopic behaviors of m suggests a coexistence of two states during the transitional periods: Representative figures of the temporal dynamics of m in both periods are presented with respect to the change of anesthetic concentration in the inset of Fig. 3D. Although the averaged behavior of the order parameter favors state 0 or 1 statistically, prolonged microscopic spatial fluctuations exist, which is a characteristics of metastability. The temporal fluctuations of up and down states in each location of the thalamocortical network are visualized in Movie S1. Taken together, these results indicate that the thalamocortical system never shuts down instantaneously in ketamine-induced loss of consciousness. Moreover, the existence of meta-stability implies that blockade and recovery of information transfer within the thalamocortical network are initiated at discrete regions fluctuating over a wide range of anesthetic concentrations, which eventually drives the whole system into one phase in a hysteretic manner.

Discussion
We defined the order parameter from experimentally obtained neuronal signals to determine the level of information blockade within the thalamocortical circuit by adopting aggregate time units and resultantly reducing the complexity of the field potential signals. Application of the aggregate time unit made it possible to split the time series into epochs containing electrophysiologically meaningful oscillations. Depending on whether each epoch showed thalamocortical resonant oscillation, we were able to determine the state vector and order parameter. Although a similar description of the system can be achieved by conventional indices such as the spectral edge frequency 95 th percentile [15,16], the discretization method established here has the advantage of tracing changes of state statistically over a biologically relevant time scale. By distinguishing each discrete epoch into an oscillatory pattern of unconsciousness or not, fluctuation between states can be discerned, and the probability of transition between states can be defined. The probability of transition between states was observed to increase in a time-locked manner to the loss and recovery of motion, indicating increased fluctuation around the behavioral transition points, and this prolonged spatial and temporal fluctuation between states implies that the transitions may not be discrete at the whole brain level. Although the anesthetic-induced transition model of Steyn-Ross et al. [5] proposed that a single macrocolumn in cerebral cortex may undergo first-order state transition with hysteresis, they also reported that the more realistic cortical model of 2-D grid of macrocolumns did not simultaneously go off to the unconscious state due to vagarious subcortical driving force [17]. Later modeling studies showed that the anesthesia-induced state transition could also occur in a nondiscrete way through Hopf bifurcation, by expanding the model to include anesthetic effects on both duration and amplitude of IPSP [6] and slow ionic currents [7]. The response curve of order parameter with respect to anesthetic concentration captured not only the path-dependent nature that the local cerebral cortical unit might have, but also the system level fluctuations during state transitions.
The investigation of the order parameter in the temporal domain revealed lag and lead between behavioral representation and the change in thalamocortical information blocking level. However, the blackout time from t LOM to t ROM demonstrated large individual variance which was mainly attributed to the variance of t ROM , raising the question that whether time is the most appropriate reference for phase transition in consciousness. Since investigation in the temporal domain might not be sufficient to determine whether the difference in the order parameter rate of change originated from pharmacokinetic or neural mechanisms, we introduced a mathematical model to infer anesthetic concentration in the brain tissue. The anesthetic concentration model was built as a threecompartment model considering the intraperitoneal injection route, where the model parameters were determined from previous pharmacokinetic studies. Basically the model parameters were estimated from the ketamine data, and the alteration of the ketamine's anesthetic effect that co-administration of xylazine could give rise was implicitly incorporated in the model parameters k a and a, which were determined to be proportional to induction time and blackout time, respectively. While the pharmacokinetics data of xylazine for mice was rarely available, co-administration of xylazine was reported to affect the action of other anesthetics in such ways of lengthening the anesthesia period in budgerigars [18] or reducing the amount of anesthetic drug used for surgical anesthesia in cats [19]. The behavioral transition points were thought as the alternative indicators of overall anesthetic effects including probably delayed recovery from ketamine anesthesia due to xylazine.
Although we adopted the behavioral transition points as proxies of anesthetic effects, the absence of simultaneously measured plasma concentration is a manifest limitation that cannot be easily evaded in estimating the anesthetic concentration. To verify the robustness of the model, its dependence on the parameter selection was tested with varying parameters of k a and a which were crucial for shaping the concentration model. As intraperitoneally injected drug is not absorbed into the circulation system directly, estimated k a , a rate constant for drug absorption from the intraperitoneal cavity, can have some variance depending on each subject's condition. While the other parameters were fixed, variation of k a up to 650% of the estimated parameter did not make qualitative changes in the order parameter curve with regard to anesthetic concentration ( Fig. S2(a) in Information S2). Variation of a, which might be raised by xylazine, resulted in wider range of variability in the order parameter curve, but the path-dependent nature was conserved even with 650% variation of a value (Fig. S2(b), Information S2). The response function of thalamocortical states with respect to anesthetic concentration demonstrated significantly lower individual variance than that with respect to time, as well as maintained its hysteretic feature regardless of variable concentration model parameters, suggesting that the anesthetic concentration is a better reference for predicting the thalamocortical system state.
The observation that a higher anesthetic concentration was required for loss of consciousness compared to the recovery period provides experimental evidence for the tendency of the brain system to resist any change in state, and supports the idea that neural activities experience neural inertia. The term ''neural inertia'' was originally introduced by Friedman et al. to describe the path-dependent behavioral responses to anesthetic dose [14]. In their work, scatter plots of behavioral responses versus anesthetic doses measured in a population demonstrated hysteresis, and Hill slopes were steeper during the induction of anesthesia compared to emergence. The smaller time derivative values of the order parameter during recovery in our study are consistent with the interpretation that the reconstitution of consciousness is more difficult than shutting it down. Our results add electrophysiological evidence that a barrier exists between the wakeful state and anesthesia-induced unconscious state even in an individual brain.
Here we have assumed that the resonant oscillation in the thalamocortical network reflects the level of information blocking, and thus the level of vigilance. The thalamocortical pathway is thought to play an important role in sensory information transfer, raising the concept of thalamic consciousness switch [20]. Corticocortical integration has also been postulated as essential for consciousness [4,21], with the experimental evidence that the corticocortical information integration level has been observed to significantly decrease during sleep and anesthesia [3,22]. Although the extent to which these mechanisms contribute for consciousness may be argued, it seems that these two mechanisms are not exclusive but rather constitutes a hierarchical structure for the emergence of consciousness. The mismatch between the blackout period and the steady state of order parameter might originate from sequential recovery of this hierarchy. We cannot say that the order parameter which quantifies the level of thalamic information blocking as the only and all-powerful index for measuring consciousness level, but we can say that the order parameter reflects the change of the primary level in the hierarchical mechanism for consciousness.
From the order parameter diagram represented as a function of anesthetic concentration, we observed the existence of metastability and path-dependency in the periods of consciousness transition induced by anesthesia. In our study, it is conceptually interpreted as the expectation value of the global resonance level among different functional units. Not restricted to this single approach, the order operator can be defined in various ways depending on the target system's characteristics and utilized independently of measurement modality and number of measurement sites. The formulation described in this work can promote the quantitative investigation of brain state transition characteristics. Unquestionably, ketamine/xylazine anesthesia is not the single route to unconsciousness, and the resonant oscillations of the thalamocortical system are not the unitary correlate of unconsciousness; however, this bottom-up approach has the advantage of enabling a gradual understanding of the ''states'' of consciousness. Figure S1 Diagram of three-compartment model. Xi denotes the amount of ketamine in the compartment Ci. The transition rates between compartments (kij) and the elimination rates from the site of drug administration (ka) and from the cerebral circulatory system (k10) are assumed to be first order. D and F mean dose of ketamine and bioavailability, respectively. The ketamine concentration at Ci is defined as the amount of ketamine at the compartment divided by the volume of the compartment.  Table S1 The parameters without asterisks are fixed ones, determined from the plasma concentration-time profiles of ketamine. Those with asterisks are adjustable ones, derived from the time profile of plasma levels of ketamine and/or the behavioral reference points i.e., t ADM , t LOM , and t ROM . The values with dagger ( { ) were expressed as means6standard deviation.

(DOCX)
Movie S1 Spatio-temporal fluctuation of states during transition periods. Upper panel shows fluctuation of the order parameter computed at the aggregate time (T) from each recording site (clockwise from top left: primary motor cortex, primary somatosensory cortex, somatosensory thalamus and motor thalamus) in the vicinity of t ADM , t LOM and t ROM . White and black colors indicate information transmission mode and information blocking mode, respectively. In the order parameter diagram at lower panel, moving red bar shows a position of corresponding aggregate time. Dashed vertical lines indicate t ADM , t LOM and t ROM , from left to right. (MP4) Information S1 (DOCX)

Author Contributions
Conceived and designed the experiments: JHC. Performed the experiments: EH. Analyzed the data: EH JHC. Contributed reagents/materials/ analysis tools: EH KH. Wrote the paper: JHC EH. Revised the manuscript: JHC SK.