Relationship of critical dynamics, functional connectivity, and states of consciousness in large-scale human brain networks

Recent modeling and empirical studies support the hypothesis that large-scale brain networks function near a critical state. Similar functional connectivity patterns derived from resting state empirical data and brain network models at criticality provide further support. However, despite the strong implication of a relationship, there has been no principled explanation of how criticality shapes the characteristic functional connectivity in large-scale brain networks. Here, we hypothesized that the network science concept of partial phase locking is the underlying mechanism of optimal functional connectivity in the resting state. We further hypothesized that the characteristic connectivity of the critical state provides a theoretical boundary to quantify how far pharmacologically or pathologically perturbed brain connectivity deviates from its critical state, which could enable the differentiation of various states of consciousness with a theory-based metric. To test the hypothesis, we used a neuroanatomically informed brain network model with the resulting source signals projected to electroencephalogram (EEG)-like sensor signals with a forward model. Phase lag entropy (PLE), a measure of phase relation diversity, was estimated and the topography of PLE was analyzed. To measure the distance from criticality, the PLE topography at a critical state was compared with those of the EEG data from baseline consciousness, iso ﬂ urane anesthesia, ketamine anesthesia, vegetative state/unresponsive wakefulness syndrome, and minimally conscious state. We demonstrate that the partial phase locking at criticality shapes the functional connectivity and asymmetric anterior-posterior PLE topography, with low (high) PLE for high (low) degree nodes. The topographical similarity and the strength of PLE differentiates various pharmacologic and pathologic states of consciousness. Moreover, this model-based EEG network analysis provides a novel metric to quantify how far a pharmacologically or pathologically perturbed brain network is away from critical state, rather than merely determining whether it is in a critical or non-critical state.

Recent modeling and empirical studies support the hypothesis that large-scale brain networks function near a critical state. Similar functional connectivity patterns derived from resting state empirical data and brain network models at criticality provide further support. However, despite the strong implication of a relationship, there has been no principled explanation of how criticality shapes the characteristic functional connectivity in large-scale brain networks. Here, we hypothesized that the network science concept of partial phase locking is the underlying mechanism of optimal functional connectivity in the resting state. We further hypothesized that the characteristic connectivity of the critical state provides a theoretical boundary to quantify how far pharmacologically or pathologically perturbed brain connectivity deviates from its critical state, which could enable the differentiation of various states of consciousness with a theory-based metric.
To test the hypothesis, we used a neuroanatomically informed brain network model with the resulting source signals projected to electroencephalogram (EEG)-like sensor signals with a forward model. Phase lag entropy (PLE), a measure of phase relation diversity, was estimated and the topography of PLE was analyzed. To measure the distance from criticality, the PLE topography at a critical state was compared with those of the EEG data from baseline consciousness, isoflurane anesthesia, ketamine anesthesia, vegetative state/unresponsive wakefulness syndrome, and minimally conscious state.
We demonstrate that the partial phase locking at criticality shapes the functional connectivity and asymmetric anterior-posterior PLE topography, with low (high) PLE for high (low) degree nodes. The topographical similarity and the strength of PLE differentiates various pharmacologic and pathologic states of consciousness. Moreover, this model-based EEG network analysis provides a novel metric to quantify how far a pharmacologically or pathologically perturbed brain network is away from critical state, rather than merely determining whether it is in a critical or non-critical state.
level (Beggs and Plenz, 2003;Friedman et al., 2012), but also at the large-scale or global network level (Eguíluz et al., 2005;Hahn et al., 2017;Haimovici et al., 2013). Until now, most studies have focused on scale-free behavior, showing power law distribution of empirically observed variables as the evidence of criticality. Recent studies have revealed that the functional connectivity of spontaneous dynamics is highly correlated to the structural connectivity when the system is at criticality, and less correlated when the system is at sub-or super-criticality Stam et al., 2016;Tagliazucchi et al., 2016). However, despite the obvious implication of a relationship between network structure, connectivity, and criticality, there has been no principled explanation of how a correlation between functional and structural networks emerges at a critical state.
Here, we hypothesized that the network science concept of partial phase locking is the underlying mechanism of the emergence of functional connectivity from a structural brain network. In a coupled oscillator model, the partially locked state emerges at an intermediate level of coupling strength, in which phases of some nodes are locked with each other while those of others drift (Ko and Ermentrout, 2008;Kuramoto, 1984). Therefore, the extent to which the phase locking shows heterogeneity in the brain network is thought to be closely associated with the global pattern of functional connectivity. Moon, Lee and colleagues identified a mathematical relationship between the phase locking and node degree of a structural network but the relationship with criticality was not explicitly tested (Moon et al., 2015). Considering the computational and empirical evidence that pharmacologically perturbed brains significantly deviate from the critical state, it may be useful to develop a metric to quantify the distance of a perturbed brain from criticality Tagliazucchi et al., 2016). Thus, we hypothesized that the distance from the critical state of a functional brain connectivity configuration of a second order phase transition provides a theory-based metric quantifying how far a pharmacologically or pathologically perturbed brain is away from critical state, which allows the comparison of various brain perturbations with the same theoretical reference.
To test the hypothesis, we used the Kuramoto model to simulate the functional connectivity in an anatomically informed human brain network model. To compare predictions from the model with empirical EEG data, we first projected the source signals of the 3-dimensional network model onto the scalp to generate EEG-like sensor signals using a forward model. By changing the control parameter, we studied the relationship between functional connectivity based on EEG and criticality. According to network science, the partially locked network at a critical state modulates both synchronization and frequency in a network heterogeneously; the higher the degree of the node, the more synchronous and slow the frequency (Moon et al., 2015). To capture these characteristic changes, we used phase lag entropy (PLE), a measure of phase relation diversity , and constructed the topographies of PLE from both the model and empirical data. The PLE topography of the model at the critical state was compared with that derived from the EEG data acquired during baseline consciousness (n ¼ 73), isoflurane anesthesia (n ¼ 29), ketamine anesthesia (n ¼ 15), vegetative state/unresponsive wakefulness syndrome (n ¼ 29), and minimally conscious state (n ¼ 16); we then measured the distance of each state from the critical state. This model-based EEG network analysis revealed a novel metric that quantifies how far a pharmacologically or pathologically perturbed brain is away from criticality, which enabled us to define different states of consciousness with a common theoretical reference point. The overall analytical scheme of this study is summarized in Fig. 1.

Criticality, partial locking, and connectivity
In this study, the critical state of the brain network was determined using two measures: (1) the large variance of order parameters (global phase synchronization) measured by the pair correlation function and (2) the large correlation between functional (EEG) and structural (anatomically informed) brain networks. The large variance of order parameter originates from metastability at a critical state (Cabral et al., 2014;Shanahan, 2010). However, the reason for a large correlation between Fig. 1. Schematic illustration of analysis 78 time series simulated by the Kuramoto model in an anatomically informed human brain network represent source signals of 78 cortical regions. The source signals were then projected onto the scalp surface using a forward model to generate EEG-like sensor signals. Empirical EEG signals were collected from three different experiments; (1) baseline (BL) vs. isoflurane anesthesia (ISO), (2) baseline (BL) vs. ketamine anesthesia (KET), (3) minimally conscious state (MCS) vs. unresponsive wakefulness syndrome (UWS). PLE was measured from both simulated and experimental signals, and topographic patterns of PLE were compared. functional and structural network has yet to be clarified. Here, we argued that the partial phase locking in network synchronization is a mechanism to shape the functional connectivity similar to the anatomical brain network at a critical state. Thus, the terminology 'partial locking' is used to reflect the underlying mechanism to explain the relationships between the functional network, structural network, and critical state, rather than just a mixed state of synchronized and incoherent connectivity. During a partial locking state, the synchronization of a node is determined by its node degree (i.e., number of connections) and local network structure; a larger degree node in the structural network has a higher synchronization with the neighboring oscillators. However, theoretically, the opposite relationship (i.e., larger node degree, lower synchrony) is also possible depending on the type of interaction function between oscillators (Ko and Ermentrout, 2008). In previous studies, we found that the relationship of larger node degree/higher synchrony holds for diverse brain networks such as human, monkey, and mouse (Moon et al., 2017). The analytic and computational approaches revealed that the relationship of larger node degree/higher synchrony/slower frequency holds for a coarse-grained spatiotemporal scale (>seconds in time, and~64 or 128 EEG channels in space) (Moon et al., 2015). In this study, we focused on large temporal and spatial scales (10-s and 21-channel EEG) that can be applied to EEG data analysis and fits well with the standard Kuramoto model. We will show how the relationship between the functional and structural networks changes as the brain networks deviate from the critical state with various perturbations such as general anesthesia and traumatic injuries.

Source signal simulation: Kuramoto model on human brain network
The Kuramoto model was implemented on a structural brain network based upon the data of (Gong et al., 2009). The study employed diffusion tensor imaging tractography to construct a macroscale anatomical network capturing the underlying common connectivity pattern of cerebral cortex among 80 healthy subjects. The 78 nodes of structural brain network were defined by 78 AAL (automated anatomical labeling (Tzourio-Mazoyer et al., 2002);) template. For each subject, cortical connection between two nodes was deemed to exist if a white matter fiber bundle existed between the two brain regions. The consistent cortical connections (a binary network) across subjects were then obtained by employing a nonparametric one-tailed sign test. Thus, the connection matrix has 78 by 78 elements; the connection between node i and j, or an element a ij of the connection matrix is 1 if two nodes are connected and 0 if disconnected.
Next, we used a Kuramoto-type model with time delay as follows (Jeong et al., 2002;Yeung and Strogatz, 1999;Zanette, 2000), where, θ j and ω j are phase (angle) and natural frequency of j th oscillator.
N is the number of oscillators (nodes), and K is coupling strength, which serves as a control parameter in the model. In the simulation, the time delay (τ i;j ) between node i and j was given proportional to the physical distance between two nodes, with propagation speed of 9 m/s. The signal of the j th oscillator (x j ) was obtained by taking the sin function of the phase θ j . The 78 signals of the model represent cortical source signals.
Frequency of node j, f j ¼ ω j =2π was randomly assigned with normal distribution (Mean AE SD: 10 Hz AE 1.5 Hz), and initial values of θs were generated from uniform random distribution ½ À π; π. We changed K from 0 to 16 with an interval of 0.1. For each K (total 161 number of Ks), 100 trials were generated with different initial values of θ and ω. For each coupling strength and for each trial, time series of 15,000 samples were simulated with sampling rate 1000 Hz. The first 5000 iterations were discarded in order to exclude a non-stationary state. The ordinary differential equations were solved by using the Runge-Kutta 2nd order method.

Sensor signal simulation via forward model
The signals generated from the Kuramoto model and structural brain network represent source activities of the brain, which are under the surface of where EEG signals are measured. In reality, the electrical potentials generated by the neural activity in the brain conduct outwards through brain tissue, the skull and finally appear at the scalp surface. In order to compare experimental EEG and the model signals, we generated surface level signals from the simulated source signals. More specifically, the 78 oscillator time series representing the source activity were converted into 21 EEG sensor level signals, which overlapped with experimental EEG channels analyzed empirically.
We used three concentric spherical head models; the three layers consist of the brain, skull, and scalp. The conductivity of the three layers was set to be 0.33, 0.0042, and 0.33 S/m, respectively (Wen and Li, 2006). The source activity was represented as a dipole moment. The coordinate of the dipole moment in the brain was determined by the region's standard coordinates and the orientation of the dipole moment was randomly assigned. We also tested the radially oriented dipoles; the results are almost the same with the random orientation of the dipole moment (Fig. S1). The forward model simulation was conducted by using the Field Trip Toolboox (Oostenveld et al., 2011).

Experimental protocol and EEG acquisition
Experiment 1: Isoflurane anesthesia The isoflurane study included 60 healthy volunteers (20-40 yrs) with body mass index less than 30. The study has been reviewed in accordance with the recommendations of the Institutional Review Boards specializing in human subject research at the University of Michigan, Ann Arbor (Protocol #HUM0071578, n ¼ 20), University of Pennsylvania (Protocol #818401, n ¼ 20), and Washington University in St. Louis (Protocol #201308073, n ¼ 20). Written informed consent in accordance with the Declaration of Helsinki was obtained from all participants. The data have been analyzed with different hypotheses and analyses Kim et al., 2018;Maier et al., 2017).
Among 20 participants at each institution, 10 of them underwent general anesthesia and 10 did not. Thus, a total of 30 participants across the three institutions were anesthetized, initially receiving propofol at increasing infusion rates over three consecutive 5-min blocks (block 1: 100 μg/kg/min, block 2: 200 μg/kg/min, block 3: 300 μg/kg/min). Loss of consciousness was measured by response to the verbal command ("Squeeze your left/right hand twice" with left/right randomized) every 30 s. Isoflurane anesthesia was then administered with air and 40% oxygen at 1.3 age-adjusted minimum alveolar concentration. Ondansetron 4 mg was administered 30-min prior to cessation of anesthesia for nausea and vomiting prophylaxis. After 3 h of isoflurane administration, the anesthetic was discontinued and responsiveness to the same verbal command was assessed every 30 s until recovery of consciousness; EEG acquired during the recovery period was not used for the current study. Considering the hypothesis, EEG data from only two states, baseline consciousness (2-min of eye-closed resting period before anesthesia) and isoflurane-induced unconsciousness (2-min of unconscious period during isoflurane administration) were used for the analysis. Because this study considers the phase relationship among oscillations, periods showing isoelectric or burst-suppression pattern were excluded from the epoch selection.
All EEG data was acquired with 32, 64, or 128-channel HydroCel nets, Net Amps 400 amplifiers (Electrical Geodesic, Inc., USA). Sampling frequency was 500 Hz and referenced to a vertex. The number of EEG channels used in the University of Michigan was 128 (n ¼ 28) or 64 (n ¼ 2), while the other two institutions used 32-channels with the same EEG Montage. We attempted to identify 21 common channels for the analysis. If there were no overlapping channels, the closest channel was chosen instead. If two electrodes were equally close, we took the average of the signals from the two (Fig. S2).
Experiment 2: Ketamine anesthesia The ketamine study included 15 healthy volunteers (20-40 yrs) with a body mass index less than 30. This study was approved by the University of Michigan Medical School Institutional Review Board, Ann Arbor, Michigan (HUM00061087), and written informed consent was obtained from all participants. All study procedures were conducted at the University of Michigan Medical School, Ann Arbor, Michigan. The data have been published with different hypotheses and analyses (Vlisides et al., 2018(Vlisides et al., , 2017. EEG data was collected during (1) 5-min eye-closed resting period before ketamine administration, (2) subanesthetic ketamine infusion (0.5 mg/kg) over 40 min with eye-closed, followed by 8 mg ondansetron, (3) break for completion of questionnaire, (4) anesthetic (1.5 mg/kg) bolus dose with eyes closed, (5) recovery period with eyes closed. Considering the hypothesis of the current study, EEG data from only two states, baseline consciousness (2-min of eye-closed resting period before ketamine administration from (1)) and ketamine-induced unconsciousness (2-min of unconscious period from epoch (4)) were used for the analysis.
The EEG data were acquired with 128-channel HydroCel nets, Net Amps 400 amplifiers (Electrical Geodesic, Inc., USA). The EEG was sampled at 500 Hz with the vertex reference. As with experiment 1, the 21 overlapping channels among 128-channels were chosen for the analysis to maximize comparison with the structural/EEG brain model.
Experiment 3: Disorders of consciousness EEG data from 80 patients suffering from disorders of consciousness due to subarachnoid hemorrhage, intracerebral hemorrhage, subdural hematoma, ischemic stroke, traumatic brain injury, meningitis or hyperglycemic brain injury were recorded on two systems. The patients were classified as either minimally conscious or in the vegetative state/unresponsive wakefulness syndrome (UWS) according to the Coma Recovery Scale Revised (Giacino et al., 2004). 17 patients (the Munich cohort) were recorded with a 64-electrode system with ring-type sintered, nonmagnetic Ag/AgCl electrodes (Easycap, Herrsching, Germany) and two 32-channel, nonmagnetic, battery-operated electroencephalographic amplifiers (BrainAmp MR, Brain Products, Gilching, Germany). The signals were recorded at 5 kHz sampling rate (BrainVision Recorder, Brain Products). 63 patients (the Burgau cohort) were recorded with a 256 channel high-density geodesic sensor net, a Net Amps 300 amplifier and Net Station 4.5. software (Electrical Geodesic Inc., Eugene, OR, USA). The signals were recorded at 1 kHz. EEG with suppression pattern gives rise to a high PLE value, which clearly results from an increased randomness of neural dynamics, not from a complex phase lead-lag relation . Therefore, we excluded data of 22 patients whose EEG shows the suppression pattern; all 22 patients were from UWS group. As with experiment 1 and 2, the 21 overlapping channels were chosen for the analysis to maximize comparison with the structural/EEG brain model (Fig. S2).

EEG preprocessing
We gathered EEG from five different states of consciousness from three experiments. Total number of data after artifact rejection was: baseline consciousness (BL, n ¼ 73), isoflurane-induced unconsciousness (ISO, n ¼ 30), ketamine-induced unconsciousness (KET, n ¼ 15), minimally conscious state (MCS, n ¼ 15), and unresponsive wakefulness syndrome (UWS, n ¼ 27). The baseline state was obtained from isoflurane (n ¼ 58) and ketamine experiments (n ¼ 15), separately. EEG from two subjects in BL state in the 1st experiment, and EEG epochs from 16 subjects in the 3rd experiment were excluded due to severe noise based on visual inspection and automatic epoch rejection.
EEG data were visually inspected and noisy periods were excluded from the epoch selection. Epochs containing amplitudes larger than 250 μV were also rejected. Data were down-sampled to 250 Hz regardless of original sampling rate, and re-referenced to the average of 21 EEG signals. For the MCS and UWS data, independent component analysis was employed such that the ocular artifacts were removed. For each state, 2min of clean data (by visual inspection) were collected, and the data subdivided into 10-s small epochs; thus, 12 epochs per each state per each subject were gathered for the PLE analysis. We applied relatively wide bandpass filtering because both pharmacologically-and pathologically-altered states are characterized by strong oscillations at relatively low frequency regime (delta~beta) and show dynamic transitions in the spectrogram. To capture the dominant oscillations of brain activities across different states, and to avoid contamination of electromyogram and electrocardiogram artifacts, we chose a frequency range, 2-25 Hz; Zero-phase bandpass filtering with FIR filter of order 1000 was conducted. Results with 2-40 Hz (broader) and 4-20 Hz (narrower) showed qualitatively similar results.

Analytic measures
Phase lag index We applied phase lag index (PLI) (Stam et al., 2007) to quantify phase-locking between two signals. Distinct from common measures of phase synchronization (Lachaux et al., 1999;Mormann et al., 2000), PLI considers only the phase lead-lag relationship of two signals. That is, its value does not depend on the real part of the cross spectrum, mitigating the effect of volume conduction in EEG recording. For large coupling strengths, consistent phase locking among signals results in a high value of PLI. On the other hand, small coupling strengths as well as independent or poorly-correlated signals result in a low PLI. The PLI between two signals can be estimated by where T is the number of time points within a given epoch, Δθ t is the phase difference between the two signals. By definition, PLI has a high (or low) value if the sign of Δθ t is consistent (or variable) within an epoch; e.g. asymmetric distribution of Δθ t results in a high value of PLI. The state or the phase of the system is indicated by global mean PLI; it indicates whether the total system is in a locked state (synchronization) or in a drift state (desynchronization). Global mean PLI is an average of PLI values from all connection pairs, i.e., for source model, a total of 3003 (¼ 78 Â 77/2) values of PLI are obtained. The global mean PLI was used as a surrogate measure of the Kuramoto order parameter < r> t ; here r is e ÀiθjðtÞ , where ψðtÞ is the average global phase. The modulus rðtÞ ¼ jzðtÞj or so-called order parameter represents the degree of synchronization, being equal to 0 when the oscillators' phases are uniformly distributed in [0, 2π) and 1 when oscillators have the same phase. The level of phase synchronization is determined by a time average of the order parameter after a transient period. N is total number of signals (nodes), and j is a signal index. Pair correlation function Critical point, or critical state, is determined by a point on the control parameter in which the divergence of the correlation length or a maximal susceptibility to an external field appears. In the Kuramoto model, a pair correlation function (C p ) measures the variance of the order parameters calculated during a time interval t, which is used as a surrogate measure of susceptibility and shows a bellshaped peak near the critical point (Yoon et al., 2015). C p in Kuramoto model is defined as where zðtÞ is defined in a rotating frame. C p measures the variance of cross-correlation coefficients among all EEG signals. Because the real part of zðtÞ is extremely vulnerable to volume conduction, C p cannot predict the critical point in the sensor signals and experimental EEGs. Thus, we identified the critical point by maximum of C p averaged across 100 trials in the source model.
Phase lag entropy Phase lag entropy (PLE) measures the diversity of phase lead-lag patterns among two signals . Unlike time-averaging measures of phase synchronization (e.g., PLI), PLE incorporates the temporal dynamics of the instantaneous phase time series into the phase synchronization analysis. More specifically, PLE extracts consecutive temporal patterns of the phase relationship, which span tens of milliseconds.
To calculate PLE between two signals, the phase difference is first symbolized in a binary fashion; the symbol s t ¼ 1 if Δθ t > 0 (first signal is phase leading the second signal), and s t ¼ 0 if Δθ t < 0 (first signal is phase lagging the second signal). Then, the vector S t representing the temporal pattern of the phase relationship is given by: where m and l represent pattern size and time lag, respectively. For instance, with m ¼ 3, eight patterns ('000', '001',' 010', '100', '011', '101', '110', and '111') can be generated. Finally, by applying the standard Shannon entropy formula to the distribution of the phase patterns, PLE is calculated, where 0 p k 1 is the probability of the kth phase lead-lag pattern, as estimated by calculating the fraction of time each pattern occurs relative to the total time length of a given epoch. In this study, we chose a parameter set m ¼ 5 and l ¼ 2 considering the time-lagged mutual information of phase lead-lag time series (Fig. S3). Testing various parameter sets did not change the PLE topography qualitatively.

Source level analysis
First, we investigated the influence of the structural nodal degree on the PLI and frequencies of 78 nodes. Specifically, we focused on the difference between hub and non-hub nodes. The hub nodes were defined as the top nine ranked nodes (node degree cutoff ¼ 13) in node degree of the structural network, and the remaining nodes were denoted as the non-hub nodes. For the PLI analysis, 78 by 78 PLI matrix was converted into a 78 by 1 PLI vector, i.e., a nodal PLI of each node was obtained by averaging all 77 PLI values. Then, PLI hubs (or PLI non-hubs ) was obtained by averaging the nodal PLI values of nine hub nodes (69 non-hub nodes). The difference between PLI hubs and PLI non-hubs was tested at three representative Ks (K ¼ 0.2, 2.4, and 10) with a Wilcoxon rank sum test, under a null hypothesis that the difference comes from a distribution symmetry of about zero; the distribution is made of 100 data points (100 simulations) for each K. The difference between the two was also compared across three representative Ks, under a null hypothesis that there are no changes in difference across two different Ks. The same statistical analysis was conducted for frequency difference. Because the variances of PLI and of frequency were different across the three representative Ks, we used a non-parametric Wilcoxon signed rank test instead of a parametric paired t-test. Next, we compared PLE topography and the structural node degree. As in PLI, 78 by 78 PLE matrix was converted into a 78 by 1 PLE vector. Then, an averaged topography of PLE was obtained by averaging the PLE topographies over 100 trials. Finally, Pearson correlation between the average PLE topography and node degree of structural network were estimated. This procedure was repeated for all K values. In all hypothesis tests, Bonferroni corrected p-values with values less than 0.05 was considered to represent significant difference across the states (*p < 0.05; **p < 0.01; ***p < 0.001).

Comparison of sensor signals and empirical data
The PLE topographies of the sensor model data and empirical data were compared with a method similar to that used for the source model analysis. In both the sensor model and the empirical data, a 21 by 21 PLE matrix was converted into a 21 by 1 PLE vector as in the source model. For the empirical data, the PLE topography per each 10-s window in each subject in each state was obtained. Then, Spearman correlation was used to evaluate the similarity between the PLE topography of the sensor model and that of the empirical data. The correlation coefficient was averaged across windows. This procedure was repeated for all K values of the sensor model and for all empirical data.
To investigate different states of consciousness, a similarity (Spearman correlation) between PLE topography of each conscious state in the empirical data and a reference PLE topography was calculated. The reference PLE topography represents a PLE vector from the sensor model at the critical point (K ¼ 2.4). Thus, high (or low) correlation value represents a state being close to (or far away from) the critical state, in terms of coupling strength, K. The similarity values of different states were compared to each other by using Wilcoxon signed rank test. We also calculated a mean PLI to evaluate the complexity of the brain connectivity for each state. The same procedure was conducted for mean PLE. The mean PLE for each subject for each state was obtained by averaging 210 (¼ 21 Â 20/2) PLE values across all 10-s epochs. The statistical test was conducted only among the same datasets (experiment), e.g., baseline consciousness (n ¼ 58) vs. isoflurane anesthesia (n ¼ 30), in the first experiment. In all hypothesis testing, Bonferroni-corrected p-values less than 0.05 were considered to represent significant difference across the states (*p < 0.05; **p < 0.01; ***p < 0.001).

Strong correlation between functional and structural brain connectivity near the critical point
In the brain network model, PLI increased and the mean frequency decreased monotonically as K increased ( Fig. 2A and B). We compared the PLIs of hub and non-hub nodes. At the intermediate level of K, the PLI of hubs is larger than that of non-hubs, meaning that hubs are more easily phase-locked than non-hub nodes, especially at intermediate values of K ( Fig. 2A, C). In this brain network model, we used PLI because it is robust with respect to the volume conduction problem of EEG. However, we also tested a local order parameter, a conventional measure for the Kuramoto model study, and obtained consistent results (Fig. S4). Fig. 2C describes the PLI differences at three representative K values (K ¼ 0.2, 2.4, and 10, respectively). The PLI differences significantly deviated from zero at K ¼ 2.4 and 10 (#: Bonferroni corrected p < 0.05). Also, the PLI difference at K ¼ 2.4 was significantly higher than that at K ¼ 0.2 and 10 (***: Bonferroni corrected p < 0.001).
Frequencies also showed disparity between hubs and non-hubs ( Fig. 2B and C). The mean frequencies of hubs and non-hubs were comparable at small and large Ks, whereas the frequency of hubs was lower than that of non-hubs at an intermediate K. This indicates slower dynamics of hub nodes compared to the non-hub nodes. Fig. 2D describes the frequency differences at three representative K values (K ¼ 0.2, 2.4, and 10, respectively). The frequency differences significantly deviated from zero at K ¼ 2.4 and 10 (#: Bonferroni corrected p < 0.05; number of hypotheses ¼ 5). Also, the frequency difference at K ¼ 2.4 was significantly higher than that at K ¼ 0.2 and 10 (***: Bonferroni corrected p < 0.001; number of hypotheses ¼ 5).
These findings are consistent with analytic results in our previous studies (Moon et al., 2017(Moon et al., , 2015. The two characteristic features of a critical state, i.e., stronger phase locking and reduced frequencies at hub nodes, can be captured with PLE, a measure of phase relation diversity. Thus, we expected that the PLE of hubs would be lower than that of other nodes near the critical point. Fig. 2E demonstrates the examples of source signals at small, intermediate (critical), and high K values. At small K (¼ 0.2), both the peripheral and hub nodes are incoherent, such that most signals drift (Fig. 2E left). At a critical point (K ¼ 2.4), hubs tend to be locked in terms of phases while non-hub nodes drift (Fig. 2E middle). When K is sufficiently high (K ¼ 10), both hub and non-hub nodes are locked and the system becomes globally synchronized (Fig. 2E, right).
We calculated the PLE of each node by averaging 77 PLE values; that is, a 78 by 78 PLE matrix was reduced to 1 by 78 PLE vector. Fig. 2F depicts the relationship between node degree and PLE. At low or high K, the correlation between PLE and node degree was relatively low (Pearson correlation coefficient, R ¼ À0.187, p ¼ 0.101 for K ¼ 0.2 and À0.358, p < 0.01 for K ¼ 10). However, at the critical point, a strong negative correlation was seen; hubs showed low PLE and peripheral nodes showed relatively higher PLE (R ¼ À0.806, p < 10 À18 ). Fig. 2G describes the change of mean PLE (averaged over all 3003 pairs) and correlation between nodal PLE and node degree as a function of K. The mean PLE decreased monotonically as K increased due to the phase-locking and frequency reduction effects. C p , a surrogate measure of network susceptibility, showed a bell-shaped curve as a function of K and was maximized at K ¼ 2.4. Here, we determined the coupling strength (K ¼ 2.4) as the critical point of the brain network model. The strong correlation between nodal PLE and node degree was seen near the critical point.

Anterior-posterior asymmetry observed at criticality in the model and EEG in resting state
The nodal PLE values of 78 regions of interest were used to construct a topographic pattern, i.e., the PLE values were represented in color on the surface of the brain (Fig. 3A). In Fig. 3A, the PLE topographies from three representative Ks are shown. For each topographic map, color was scaled such that minimum (mean -2SD) and maximum values (mean þ 2SD) were determined. Among the three different Ks, a notable topographic pattern was observed at criticality. At the critical point, the topography of PLE showed an anterior-posterior asymmetry; PLE is higher in frontal regions and lower in posterior regions (Fig. 3A, middle). This is because of the spatial distribution of node degree of the human brain network (Gong et al., 2009). The human brain network has many strong hub nodes in posterior regions and many peripheral nodes in anterior regions (Fig. 3B) so that the strong negative correlation between node degree and PLE at the critical point ( Fig. 2E and F) leads to the anterior-posterior asymmetry.
We further investigated whether the sensor model and experimental EEG can exhibit anterior-posterior PLE asymmetry. Fig. 3C shows PLE topography of the sensor signals from K ¼ 0.2, 2.4, and 10, respectively. The asymmetric anterior-posterior PLE was clearly seen at K ¼ 2.4, but not at K ¼ 0.2 or 10, as in the source model (Fig. 3A).
Experimental EEG recorded during the baseline conscious state also exhibited the asymmetric anterior-posterior PLE topography. Mean PLE topography is shown in Fig. 3D (n ¼ 73). We quantified the similarity between the sensor model and experimental EEG. The Spearman correlation was calculated between PLE of EEG data and that of the sensor models from different Ks. Fig. 3E shows that the Spearman correlation is higher near the critical point than low or high K values. Fig. S5 describes the same analysis result based on PLI topography. The Spearman correlation between PLI topography of EEG data and that of the sensor model was maximal near the critical point as well, but the Spearman correlation values were substantially lower than that of PLE (Maximum Fig. 2. A strong negative correlation between PLE and node degree is observed near the critical point. (A) As K increases, both PLI of hubs and non-hubs increase. The PLI of hubs increases faster than that of non-hubs, resulting in a disparity at intermediate K. (B) As K increases, the frequency (f) of both hubs and non-hubs decreases. The f of hubs decreases faster than that of non-hubs, resulting in a disparity at intermediate K. (C) The disparity in PLI is significantly larger at K ¼ 2.4 than K ¼ 0.2 and K ¼ 10. (D) The disparity in f is significantly larger at K ¼ 2.4 than K ¼ 0.2 and K ¼ 10. (E) Representative signals of hub and non-hub nodes at three different Ks. At K ¼ 2.4, signals of low degree regions change their phase lead and lag relationship frequently, whereas signals of high degree regions are locked with each other. (F) Correlation between PLE and node degree at three different Ks. A strong correlation is seen at K ¼ 2.4, due to the disparity between hub and nonhub nodes in PLI and frequency near the critical point (B-C). (G) Mean PLE, C p , and absolute value of correlation between PLE and node degree is shown as a function of K. When C p is maximized, i.e., near the critical point, the correlation between PLE and node degree is maximized. correlation ¼ 0.240 and 0.541 for PLI and PLE, respectively). PLE better reflects the disparity of hub-parietal activities at the critical point than does PLI. The topoplots of PLE and PLI with absolute scale are depicted in Fig. S6.

Disrupted partial phase locking in pharmacologically and pathologically perturbed brains
We tested how pharmacological and pathological perturbations of the human brain network alter the topographical pattern of PLE, with comparison to the model brain network. The topographic patterns of PLE for the five states of consciousness (BL, ISO, KET, UWS, and MCS) are presented in Fig. 4 (A). The topographical similarity curve was generated by comparing the topographical pattern of PLE from empirical EEG and all topographical patterns generated from the model brain network as changing the coupling strength K (from 0 to 16) (Fig. 4B); this procedure was performed for each epoch in each subject in each state. The BL states of consciousness in the two anesthetic experiments had a bell-shape curve with the maximum topographical similarity with the topographical patterns around K ¼ 2.4 (which is considered to be the critical point in the model brain network). By contrast, the two anesthetics showed similarity curves that were distinct from those of the baseline states. In both ISO and KET states, the overall topographical similarity was lower, especially without an obvious maximum peak in the similarity curve. In the topographical pattern, the frontal PLEs were very low (Fig. 4A), which is distinct from the baseline conscious state. Of note, the topographical pattern of MCS and UWS patients had a bell-shape curve, although both similarity levels were lower. The topographical similarity of both the UWS and MCS patients were not significantly different from each other. The different similarity curves between the two anesthetics, UWS, and MCS may imply distinctive effects of pharmacological and pathological perturbations on the topographical pattern of PLE. For the UWS and MCS patients, a few primary structures of the topographical pattern (for instance, the asymmetric anterior-posterior PLE) are preserved, but the secondary structures are altered, whereas the two anesthetics perturb the overall structure. Fig. 4C presents a 2-dimensional parameter space that consists of topographic similarity of PLE referenced to the model brain network at a critical point and mean PLE. The topographical similarity of PLE distinguished the BL state from all other perturbed networks. By contrast, the mean PLE differentiated the BL and MCS from two anesthetized states and UWS, respectively. In other words, the topographic similarity distinguishes the normal and the perturbed brains, whereas the mean PLE distinguishes conscious (baseline and MCS) and unconscious states (anesthetized and UWS). Both ISO and KET states showed a reduction of similarity and of mean PLE, compared to the baseline (Spearman correlation: ###p < 0.001 for isoflurane and ketamine; for mean PLE: ***p < 0.001 for isoflurane, *p < 0.05 for ketamine). The mean PLE of  In (A, C, and D), color bar presents a relative scale (; the mean AE2SD of 21channels for each K) was applied.
MCS was significantly higher than that of UWS (**: p < 0.01). However, the topographical similarity of the UWS and MCS patients were not significantly different from each other (Bonferroni unadjusted pvalue ¼ 0.834). The p-values were Bonferroni corrected across the number of hypotheses (¼ 6).

Scale dependent behavior of PLE topography
We varied the window size for the PLE calculations from 200-ms to 15-s and assessed whether there is a scale dependence of PLE topography. In this analysis, Spearman correlation between nodal PLE of the sensor model at critical point and that of the empirical data was calculated per each window and then averaged across windows.
The correlation between the reference PLE vector and nodal PLE of empirical data at the baseline states increased as a function of window size (Fig. 5). An abrupt increase of the correlation was observed in 200ms to 5-s. However, in the four altered states of consciousness, the correlation values stayed low and did not increase abruptly as in the baseline states. The baselines demonstrate more significant correlations for all window sizes and apparent scale dependence of the correlation that differentiates the baseline state from the four altered states of consciousness. A smaller (or larger) correlation implies that the PLE patterns are different from (or similar to) the reference PLE pattern. The results indicate that there is a window size that is able to associate the PLE topography with the structural connectivity in the baseline state. By contrast, the PLE topographies of the perturbed brains are not constrained by the structural connectivity at any temporal scale.

Brain networks at criticality shape functional brain connectivity
Modeling and empirical studies have provided evidence that the resting-state brain is characterized by criticality. Power law distribution, large variability, and slowing dynamics of neurophysiological data support this hypothesis (Beggs and Plenz, 2003;Deco et al., 2017;Eguíluz et al., 2005;Miller et al., 2009;Tagliazucchi et al., 2012). In particular, the maximal correlation between the functional connectivity of the empirical data and the brain network model at criticality is direct evidence suggesting that the resting-state brain functions at criticality (Deco and Jirsa, 2012). fMRI studies demonstrated that the correlation is state specific, with higher values in the conscious state and lower values in the unconscious state (Tagliazucchi et al., 2016). Although these modeling and empirical studies suggest a relationship between network connectivity, state, and criticality, the nature of this relationship had not been identified. In this study, we provide evidence that partial phase locking is the network mechanism that defines functional connectivity patterns within the scaffold of structural connectivity. Furthermore, we demonstrate for the first time that the resting-state functional connectivity of the EEG, a modality distinct from fMRI, also has a large correlation with the model data at a critical state. Because of the limitations of localizing source-level signals, EEG has not been used to study the relationship between functional and structural connectivity in the brain. In this study, by using the forward model, we projected the source signals onto the scalp to generate EEG-like sensor signals, such that direct investigation of the structure-function relationship with EEG was possible. We then compared directly the PLE topographies of both the EEG and the simulation of the sensor level, and showed that the characteristic connectivity pattern, asymmetric anterior-posterior PLE, during the conscious state represent statistical significance in mean PLE and Spearman correlation, respectively. The statistical test was performed within the same experimental dataset; p-values were adjusted by Bonferroni correction (* or #: p < 0.05, ** or ##: p < 0.01, and *** or ###: p < 0.001). was only observed when the model dynamics were near the critical point. We also demonstrated that pharmacologically or pathologically perturbed states (ketamine anesthesia, isoflurane anesthesia, UWS, and MCS) were separated by the topographic similarity of PLE and the mean PLE, which yielded distinct information for a perturbed brain network. The topographic similarity distinguishes the normal and the perturbed brains, whereas the mean PLE distinguishes conscious and unconscious states, enabling us to distinguish among diverse brain states with EEG. The model-based EEG connectivity at criticality can potentially be used as a theoretical reference point to define a perturbed brain state, which might be better than comparison to a baseline EEG that cannot always be guaranteed as the optimal brain state. The topographic similarity of PLE across subjects in each state (Baseline, ISO, KET, UWS, and MCS) also supports our hypothesis that the brain network structure at criticality shapes the PLE topography (Fig. S7). BL has higher topographic similarities of PLE across subjects within and between the baselines of two anesthesia experiments (Spearman correlation: 0.601 AE 0.185 and 0.571 AE 0.179 within 1st and 2nd experiment, respectively; 0.557 AE 0.197 between 1st and 2nd experiments). By contrast, the perturbed brain states have low or no topographic similarities of PLE across subjects (Spearman correlations, ISO: 0.356 AE 0.266, KET: 0.395 AE 0.349, MCS: 0.054 AE 0.364, UWS: 0.188 AE 0.265). The pharmacologic and pathologic perturbations significantly disrupt a common topographic structure of PLE among the subjects.

Scale dependency of the similarity of functional-structural connectivity
It is noteworthy that there is a discrepancy among fMRI studies regarding the similarity of functional-structural connectivity during conscious and anesthetized states. Several studies reported that functional connectivity of anesthetic-induced unconscious state is closer to the structural connectivity than that of conscious resting state (Barttfeld et al., 2015;Ma et al., 2017;Mashour, 2018;Uhrig et al., 2018). On the contrary (Tagliazucchi et al., 2016), demonstrated a larger similarity of functional and structural connectivity in consciousness, with the similarity diminishing in an anesthetized state. (Barttfeld et al., 2015;Uhrig et al., 2018) and (Ma et al., 2017) focused on the dynamics of functional connectivity of fMRI in a relatively small window size (<¼ 3 min), suggesting that a diverse repertoire in the resting state is a signature of consciousness to be differentiated from anesthetic-induced unconsciousness. On the other hand (Tagliazucchi et al., 2016), highlighted the similarity of functional and structural connectivity that requires a whole fMRI dataset without windowing to contain the characteristic brain dynamics, i.e., slowing frequency and long range temporal correlation, at the critical state. The discrepancy among the previous fMRI data may result from the different time scales the studies were focused on: one in a relatively short time window for dynamic repertoires and the other in a large time window for comparison of functional-structural connectivity.
In our EEG study, the similarity between PLE of empirical data and of the model was pronounced in large windows, but diminished in small windows (<5 s). Importantly, the scale-dependency appears only in the baseline conscious state, not in the altered states of consciousness. This result implies that the conscious brain is characterized by diverse repertoires of global functional connectivity on a short time scale. When combining all repertoires of small windows in a large window, the functional connectivity reflects the constraints of the structural connectivity, showing a large correlation of functional and structural connectivity. Pharmacological or pathological perturbation of brain networks may reduce the diverse repertoires of functional connectivity, with the few patterns that remain reconfiguring the PLE.

Partial phase locking as a mechanism of flexible brain connectivity
At an intermediate level of coupling strength in our model, the system is in between order (synchronization) and disorder (desynchronization). The critical state, which is characterized by a partially locked state in a coupled oscillator model, consists of phase-locked subpopulations and phase-drifting subpopulations (Ko and Ermentrout, 2008;Kuramoto, 1984). For some nodes (i.e., the phase locked subpopulation), the coupling strength is strong enough to promote synchronization, but for the rest of the nodes, it is not strong enough and the phases incoherently drift with reference to each other. Studies of coupled oscillator models have shown that the partially locked state can emerge from a heterogeneity of intrinsic frequencies or the heterogeneity of node degree (Acebr on et al., 2005;Ko and Ermentrout, 2008;Strogatz, 2000). In our study, the partially locked state resulted from a disparity in node degree, giving rise to hub nodes in a synchronized state and non-hub nodes in an incoherent state. In addition, a transition from drift (incoherent) to locked (synchronous) state accompanies the slowing-down of frequencies ( Fig. 2A-C) (Moon et al., 2015). That is, when nodes become synchronous, they oscillate with slower frequencies compared to their intrinsic frequencies. Therefore, hubs have a higher probability of being in a locked state with slower frequencies, whereas non-hub nodes are likely in a drift state with relatively faster frequencies. PLE is inherently sensitive to changes in both phase-locking and frequencies of oscillators , thereby showing a pronounced negative correlation between nodal PLE and node degree at the intermediate coupling strength. We defined the intermediate coupling strength as the critical point. In the human brain network, the partially locked state maximizes the influence of the structural connectivity on the functional connectivity, and the broad ranges of functional connectivity, frequency, and power at the critical point induce susceptibility to perturbation in the network. If a system deviates from the critical point, it begins to be biased toward an incoherent or synchronized state, with loss of variability and loss of susceptibility.

Criticality and altered states of consciousness
If the conscious brain operates near the critical point, what happens in pharmacologically or pathologically perturbed brains? What about other altered states of consciousness such as sleep or epileptic seizures? Hudetz et al. used a spin-glass model to test whether anesthesia or epileptic seizure would move the brain away from the critical point and reduce the diversity of connectivity patterns . A recent EEG study reported that the conscious state is characterized by a high dynamic complexity, implying a balanced state between order and disorder . Disruption of scale-free organization of neural activity has been observed in various altered states of consciousness with different sources of data from different species. In nonhuman primates, both propofol-and ketamine/medetomidine-induced loss of consciousness were associated with stabilization of cortical dynamics (Solovey et al., 2015). During epileptic seizures, neuronal activity patterns deviate from power law distribution (Meisel et al., 2012). During emergence from pentobarbital anesthesia in mice, scale-invariant spatiotemporal patterns of neural activity gradually emerged based on voltage imaging (Scott et al., 2014). Tagliazucchi et al. revealed that the correlation between structural connectivity and functional connectivity of fMRI data was significantly reduced in propofol-induced unconsciousness, and argued that anesthetic-induced unconsciousness is associated with a departure from critical dynamics (Tagliazucchi et al., 2016). In our study, the unconscious brain, either induced by anesthesia or a disorder of consciousness, also showed distinct PLE topography compared to that of the conscious brain. However, the different PLE patterns cannot completely rule out the possibility that the unconscious brain might still operate near the critical point. For instance, anesthesia may reconfigure the brain network (Lee et al., 2013;Schroter et al., 2012), and, as a result, a different connectivity pattern may arise that is still governed by critical dynamics despite a perturbed condition. Additionally, there were empirical studies that observed a preservation of scale-free organization of EEG-derived functional connectivity during propofol-induced unconsciousness (Hahn et al., 2017;Lee et al., 2010). Furthermore, neural avalanche distributions from human local field potential recordings followed a power law during rapid and non-rapid eye movement sleep as well as wakefulness (Priesemann et al., 2013). Liu et al. found that the scale-free distributions of node size and node degree in fMRI were preserved across wakefulness, propofol sedation, and recovery. However, the scale-free distribution was notably absent in UWS patients (Liu et al., 2014). The results imply that general anesthesia does not seem to be a complete network failure but rather that the brain undergoes an adaptive reconfiguration to maintain an optimized state of global brain network organization.

Application to altered states of consciousness
PLE measures the entropy of phase lead-lag patterns between two EEG signals . Higher PLE reflects more complex interactions between brain regions, which correspond to expanded repertoires in inter-areal communication. In contrast, the topographical similarity of PLE measures how close the global communication structure derived from empirical data is to that of the theoretically estimated critical state. In the 2-dimensional state space (Fig. 4C), we observed that the mean PLE separates the baseline conscious states and MCS from anesthetized states and UWS, and the topographic similarity of PLE differentiates the baseline conscious states from all other perturbed brains (anesthetized states, MCS, and UWS). The baseline showed a larger mean PLE and higher topological similarity compared to the theoretically optimal topographic pattern. Two anesthetized states and UWS had smaller mean PLE and lower topological similarity. Notably, MCS has a lower topological similarity but a higher mean PLE. The results suggest that the overall communication complexity of the MCS brain is at a similar level with the baseline, but the communication structure is inefficient and far from the critical state, thus reducing sensitivity to external stimuli. These PLE based measurements are consistent with the perturbational complexity index (PCI), which quantifies the complexity (algorithmic compressibility) of the EEG response to a direct cortical perturbation with transcranial magnetic stimulation (Casali et al., 2013;Sarasso et al., 2015). PCI has been proposed as the most promising index for quantifying the level of consciousness and shows high accuracy in detecting consciousness in a large population of subjects. However, PCI requires external magnetic stimulation to the brain and is a data driven index; thus, it does not have a theoretical boundary. By definition, PCI as well as PLE measure the randomness of a system but this does not allow the estimation of the criticality of a system, i.e., how far the given system deviates from its optimal state. As such, the 2-dimensional state space described in this study could be useful to define various states of normal and perturbed brain networks with a theoretical reference point.

Limitations
The study has notable limitations. First, it was a population study and thus not appropriate for individual assessment that would be more clinically relevant. However, as shown in Fig. S8, we conducted a regression analysis assessing the similarity of PLE, mean PLE, and the CRS-R (Coma Recovery Scale-Revised) scores of UWS/MCS patients (Table S1). The topographic similarity of PLE did not correlate with the CRS-R scores (R ¼ 0.082, p > 0.05), whereas the mean PLE (less affected by the network structure) was significantly correlated (R ¼ 0.461, p ¼ 0.002). Since the PLE topography of the model network was constructed based on the brain network structure of a healthy young subject population (18-31yrs), it was not appropriate to compare topographical similarities with heterogeneous individuals in the UWS/MCS group with a broad range of ages (53.27 AE 15.83 yrs; 19-90 yrs), different etiologies, and varying lesions in the brain networks. We predict that matching individual brain network structure with the same individual's functional/ neurophysiological data will improve the reliability of regression analysis.
Second, we did not examine the temporal evolution of the topographic similarity of PLE, which could enable us to test its potential as a metric. We added two examples of the temporal evolutions of both the topographic similarity of PLE and the mean PLE during general anesthesia (Fig. S9). They correlated well with state transitions during the loss and recovery of consciousness. Further work is required to associate both indices with the anesthetic depth of an individual.
Third, in order to compare the multi-institutional EEG data, we used only 21 common EEG channels, which limits spatial resolution. In the forward model, the projection of source signals (N ¼ 78) to sensor signals (N ¼ 21) may constrain the performance.
Finally, when applying time delay (τ i;j ) to the Kuramoto model, we used Euclidean distances between two areas instead of the actual length of the streamlines. In the additional analysis, we applied a fixed time delay, which did not change the results qualitatively as long as the time delay is less than a quarter cycle (<25 ms) of the natural frequency.

Conclusion
The modeling results and empirical EEG data from normal and perturbed brain networks demonstrated that the optimal functional connectivity emerges near the critical point. The model brain network revealed that partial phase locking is a mechanism of optimal functional brain connectivity. The functional connectivity emerging from the partially locked state at criticality can be used as a theoretical reference to define a perturbed brain state, which allows the quantification of how far a perturbed brain network is from its optimal functional connectivity. Our approach based on EEG and modeling may provide a novel method to monitor brain state transitions and further understand the relationship of criticality, connectivity, and states of consciousness.

Conflicts of interest
We declare no conflict of interest for all authors.