Criticality as a Determinant of Integrated Information Φ in Human Brain Networks

Integrated information theory (IIT) describes consciousness as information integrated across highly differentiated but irreducible constituent parts in a system. However, in a complex dynamic system such as the brain, the optimal conditions for large integrated information systems have not been elucidated. In this study, we hypothesized that network criticality, a balanced state between a large variation in functional network configuration and a large constraint on structural network configuration, may be the basis of the emergence of a large Φ¯, a surrogate of integrated information. We also hypothesized that as consciousness diminishes, the brain loses network criticality and Φ¯ decreases. We tested these hypotheses with a large-scale brain network model and high-density electroencephalography (EEG) acquired during various levels of human consciousness under general anesthesia. In the modeling study, maximal criticality coincided with maximal Φ¯. The EEG study demonstrated an explicit relationship between Φ¯, criticality, and level of consciousness. The conscious resting state showed the largest Φ¯ and criticality, whereas the balance between variation and constraint in the brain network broke down as the response rate dwindled. The results suggest network criticality as a necessary condition of a large Φ¯ in the human brain.


Introduction
Integrated information theory (IIT) proposes that consciousness equates with integrated information in a system and the integrated information is maximized when integration and differentiation of the systems' components are balanced [1][2][3][4][5][6]. However, in a dynamically networked system such as the brain, the optimal conditions under which the integrated information is maximized has not been elucidated. Since the integrated information in the brain network is significantly associated with functional connectivity between brain regions, it is important to identify network conditions that can facilitate the development of large amounts of integrated information. In this study, we hypothesized that network criticality, a balanced state between a large variation of functional network configurations and a large constraint of structural network configurations, may be the basis of the high Φ in conscious brains.
Criticality was originally introduced for studying phase transition in physics, which was simply defined as a balanced state between order and disorder in the activities of the elements that make up a system [7,8]. This property has been observed broadly in physical and non-physical systems and has been suggested as an optimal state for information storage, transmission, and integration with high susceptibility to external stimuli [9][10][11][12][13]. In particular, several computational modeling and empirical studies suggest that the brain dynamics associated with consciousness reside near A simple coupled oscillator model (Kuramoto model) was implemented on an anatomically informed human brain network structure constructed from diffusion tensor imaging (DTI). The Φ was calculated while modulating the level of criticality with a control parameter (coupling strength). The level of criticality was defined using PCF, a susceptibility measure, of the simulated brain activity. (B) Empirical study: 64-channel electroencephalography (EEG) data derived from seven healthy volunteers were recorded while gradually increasing sevoflurane concentration from 0.4% to 0.6% to 0.8%, then decreasing it from 0.8% to 0.6% to 0.4%. The changes of PCF and Φ were compared with the response rate to verbal commands.
To identify a relationship between criticality and Φ, we analyzed computationally a large-scale human brain network model adjusting criticality with a control parameter ( Figure 1A). The criticality A simple coupled oscillator model (Kuramoto model) was implemented on an anatomically informed human brain network structure constructed from diffusion tensor imaging (DTI). The Φ was calculated while modulating the level of criticality with a control parameter (coupling strength). The level of criticality was defined using PCF, a susceptibility measure, of the simulated brain activity. (B) Empirical study: 64-channel electroencephalography (EEG) data derived from seven healthy volunteers were recorded while gradually increasing sevoflurane concentration from 0.4% to 0.6% to 0.8%, then decreasing it from 0.8% to 0.6% to 0.4%. The changes of PCF and Φ were compared with the response rate to verbal commands. To identify a relationship between criticality and Φ, we analyzed computationally a large-scale human brain network model adjusting criticality with a control parameter ( Figure 1A). The criticality was defined using the pair correlation function (PCF), a surrogate measure of susceptibility [24,25]. The relative change of Φ across states was estimated with Φ, a surrogate measure of Φ that we developed for high-density electroencephalography (EEG; Figure 1B). Empirically, we modulated the level of human consciousness in a stepwise manner with an anesthetic and calculated both Φ and PCF from continuous EEG data. From the modeling study and EEG analysis, we were able to quantitatively study the relationship between Φ, criticality, and consciousness, and suggested network criticality as a necessary condition for a high Φ in human brain networks.

Methods
Ethics Statement: This study was conducted at the University of Michigan Medical School and approved by the Institutional Review Board (HUM00061087); after careful discussion, written informed consent was obtained from all participants. Human brain network model: Many recent studies have successfully applied Kuramoto/Stuart-Landau models to the brain in order to understand the organizational principles of multiscale brain function, surrogates of information flow, and complex dynamics at the whole brain network level [16,26,27]. Similarly, we hypothesized that the application of simple oscillatory models, which can modulate the criticality with a control parameter, to an anatomically informed brain network structure could inform the relationship between brain network criticality and integrated information.
We used a large-scale brain network model that implements a coupled simple oscillator model on the scaffold of an anatomically informed human brain network structure. The human brain network consists of 78 parcels of the cerebral cortex constructed from diffusion tensor imaging (DTI) of 80 young adults [28].
Here S is the coupling strength between oscillators and A jk denotes the anatomical connections between oscillator j and k, yielding 1 if a connection exists and 0 otherwise. τ jk is the time delay between node j and k. θ j (t) is the phase of oscillator j at time t. ω j is the intrinsic frequency of oscillator j. The Kuramoto model is the canonical model for coupled oscillators as a first-order approximation of more complex couple oscillatory systems. In general, at a sufficient coupling strength, a system of near-identical coupled limit-cycle oscillators can be approximated by this general phase model [29]. In other words, simulation results with the Kuramoto model hold for other complex network models when simplified with first-order approximation or mean-field approximation. Considering that EEG reflects superficial activity, the EEG network is a good case for the application of the Kuramoto model.
Simulation procedures: All parameters for the models were set to simulate alpha oscillations in the brain. Alpha oscillation models have successfully explained empirically observed brain network behaviors such as functional connectivity, traveling waves, and network state transitions based on electroencephalography (EEG) and magnetoencephalography (MEG) [16,21,[30][31][32][33][34][35][36][37]. Thus, alpha oscillations were analyzed to understand the behaviors of Φ and the criticality at the brain network level. The natural frequencies of the oscillators in our simulation were given as a Gaussian distribution around 10 Hz with a standard deviation of 0.5 Hz. Time delay was set proportional to the physical distances between nodes with a propagation speed of 8.6 m/s [38]. The coupling strength between the oscillators was continuously changed from 0 to 18 with an increase of 0.2 within a single configuration. For each coupling strength, we calculated criticality and integrated information. One hundred configurations were simulated with randomly selected initial frequencies and phases. The results were averaged over all the configurations.
Experimental procedures: We conducted a secondary analysis of high-density EEG data from a study of sevoflurane-induced unconsciousness in humans; the details of the experiment can be found in the supplementary text of this article and the previously published Blain-Moraes et al. [39]. The current study tested different hypotheses and, unlike the original study, included computational analyses of brain network models.
Sevoflurane dataset: 64-channel EEG were recorded from seven healthy volunteers as sevoflurane concentrations in high-flow oxygen (8 L/min) were gradually increased from 0.4% to 0.6% to 0.8% (the average range at which unconsciousness was induced) or beyond, then decreased from 0.8% to 0.6% to 0.4%. The EEG was recorded with eyes closed. The loss and recovery of consciousness were defined as the loss and recovery of response to the verbal command 'squeeze your left (or right) hand twice,' on a recorded loop every 30 s, with right/left hand commands randomized.
Due to the heavy computational load of integrated information, we did not use the whole EEG recording during the experiment. To reduce computation time, we chose the first 12-second-long EEG epoch of every minute, which was enough to represent the continuous changes in integrated information and criticality. Considering the important role of the alpha band (8)(9)(10)(11)(12) in the global information integration [16,21,[30][31][32][33][34][35][36][37], we also focused on the alpha frequency band (8)(9)(10)(11)(12) in the EEG analysis. The average reference was used for re-referencing and the windowed sinc-FIR filter (in the MATLAB toolbox from EEGLAB) was used to avoid a possible shifting of the signal.
Criticality: Criticality, a boundary state between order and disorder, has long been proposed to play an important role in neural dynamics and brain function. Empirical evidence supports the hypothesis that the brain operates at or near the critical point, not only at the neuronal network level [22,40,41], but also at the large-scale or global network level [9,12,[42][43][44]. Until now, most studies have focused on scale-free behavior, showing power law distribution of empirically observed variables. It has also been recently proposed that high correlation between functional and structural brain networks [14,16,23,26,45] and a large pair correlation function (PCF) [24,46] is evidence of criticality. In both the brain network model and empirical EEG data, we estimated criticality with PCF, which is the variance of global phase synchronization and is equivalent to susceptibility in statistical physics: where Re[z(t)] is the real part of the z(t) in Equation (3).
where Ψ is the order parameter phase. The absolute value r = |z| represents the degree of synchronization. The r is equal to zero when the phases of nodes are uniformly distributed and one when all the nodes have the same phase. Calculation of Φ: IIT defines integrated information (Φ) as the effective information (ϕ) of the minimum information partition (MIP) in a system [1][2][3]5,47]. The MIP is also defined as the partition having minimum effective information among all possible partitions.
MIP (x) = arg min ϕ (X; x, P) , where X is the system, x is a state, and P is a partition P = M 1 , . . . , M r . Identifying the MIP requires searching all possible partitions and comparing their effective information ϕ. This is the most time-consuming and computationally demanding process in the application to high-density EEG. Furthermore, considering the fact that EEG data recorded during anesthetic state transitions are non-Gaussian and continuous time series, we used the Φ AR as a measure of integrated information [48]. Φ AR is a measure of integrated information for systems with a non-Gaussian distribution of time series. The covariance of the time series of the empirical integrated information is substituted with the prediction error of the linear regression of time series.
where τ is the time delay between the past and present states, and M 1,2 are the bi-partitioned subsets. E is the prediction error of linear regression, The computed Φ AR values from original signals were compared with the Φ AR values of surrogate data sets, and non-significant Φ AR values were set as zero.
Notably, Φ AR is based on IIT 2.0 that is different from the latest version IIT 3.0. Over the last decade, the improvement of unrealistic computation time was an important issue because of the need to search an enormous number of partitions to identify the MIP. In our previous study, we proposed a surrogate measure Φ to circumvent the explosive computational time of Φ for high-density EEGs. Φ estimates the relative change of Φ across states by considering the average feature of many small sample units rather than trying to identify the MIP and its effective information for all EEG channels. A sample unit consists of a small number of EEG channels randomly selected from 64 channels and the total number of sample units is taken as large enough to represent the behavior of the entire high-density EEG montage. In this study, we limited our interest only to the relative changes of Φ values across states, rather than attempting to calculate the exact Φ value for each brain state, which would be impossible to measure using the superficial and spatially imprecise brain activity recorded by EEGs. We also verified our result by comparing it with the Φ calculated by phi toolbox, provided in Kitazono et al. [49]. In this toolbox, the Queyranne's submodular optimization algorithm was applied to reduce the number of partitions (Supplementary Figure S1).
For each sample unit, we were able to calculate the MIP and its effective information-that is, the Φ of the sample unit. For instance, in this study, we selected eight random EEG channels for a sample unit and acquired 200 sample units that were randomly selected from the baseline states. The same 200 sample units determined in the baseline were then compared across EEG windows to investigate the increase or decrease of Φ values. Since the number of possible bipartitions of eight channels is 127, calculating Φ for all EEG windows of seven subjects during state transitions is possible within a relatively short computational time.
The average Φ is defined as follows.
where n is the number of EEG channels for a sample unit and k is the number of sample units of the n EEG channels. Φ i (n) measures the effective information of MIP for the n EEG channels, by definition, the integrated information of the sample unit. Here, we chose n = 8 and k = 200 following our previous study [50]. Φ surr (n) is the spurious Φ i (n) estimated from randomly shuffled EEG data sets. Subtracting the randomness, Φ reflects the average integrated information of 200 sample units taken from the high-density EEG data that exceeds the spurious information integrated from the surrogate data. Since Φ estimates the relative change of Φ, it is appropriate for our purposes to detect the maximum Φ to compare with the maximum criticality. Functional connectivity: Phase lag index (PLI), a measure of phase locking between two signals, was used to define the functional connectivity in the EEG network [51]. We chose a Hilbert transform to extract the instantaneous phase of the EEG from each channel and calculate the phase difference ∆θ ij (t) between channels i and j, where ∆θ ij (t) = θ i (t) − θ j (t), t = 1,2, . . . ,n, and n is the number of samples within one epoch. PLI ij between two nodes i and j is then calculated using Equation (7): where, the sign() function yields: 1 if ∆θ ij (t) > 0; 0 if ∆θ ij (t) = 0; and −1 if ∆θ ij (t) < 0. The mean < > is taken over all t = 1,2, . . . ,n. If the instantaneous phase of one signal is consistently ahead of the other signal, the phases are considered locked and PLI ij ≈ 1. However, if the signals randomly alternate between a phase lead and phase lag relationship, there is no phase locking and PLI ij ≈ 0. Surrogate data: To control for spurious connectivity of EEG, 20 surrogate data sets were generated with a random shuffling method, in which a time point is randomly chosen in each EEG channel; the EEG epochs are then shuffled before and after the time point. The shuffled data have the same amplitude distribution and power spectrum of the original EEG, but there are disruptions of the original connectivity between two EEG signals.
Network construction: We expected that different EEG frequency bands and different states would have different levels of spurious connectivity [52]. Thus, after subtracting the median PLI of 20 surrogate data sets, if the remaining PLI was larger than 0.1 then the connectivity of two EEG signals was set as 1; otherwise, it was set as 0. The threshold (0.1) was chosen to avoid isolated nodes in the EEG network in the baseline states. The node degree of an EEG channel was defined as the number of functional links with other channels in the network.
EEG simulation: To test if the PLI network of EEG in conscious states is similar to the structural brain network in a critical state, we simulated 78 source signals in the structural brain network in a critical state and projected it into the 64 sensor signals on the scalp. We could then directly compare the PLI networks of EEG and the PLI network of the simulated EEG (the 64 sensor signals on the scalp). 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 is measured. In reality, the electrical potentials generated by the neural activity in the brain conduct outwards through the brain tissue and the skull and finally appear at the scalp surface where the EEG signal is measured. In order to compare experimental EEG and the model signals, we generated surface level signals from the simulated source signals. 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 [53]. 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. The forward model simulation was conducted by using the Field Trip Toolboox [54].
Statistical Analysis: We performed a one-way ANOVA ("anova1.m", MATLAB toolbox) with Tukey-Kramer correction ("multcompare.m" with alpha = 0.05 and ctype = "tukey-kramer" in MATLAB) for the comparison of the PCF and Φ among different states. The statistical tests were carried out for the modeling and empirical analysis separately. The adjusted p-values of 0.05 or lower (* p < 0.05, ** p < 0.01, and *** p < 0.001) were considered to be statistically significant (Tables S1 and S2).

Results Model Study: Correlation Between Criticality and Φ
We simulated various brain network dynamics by adjusting the control parameter (i.e., coupling strength). First, we determined the critical state of the model by identifying the coupling strength that yielded a maximum criticality. The criticality of the system was defined by PCF, which reflects the susceptibility of the brain network dynamics to perturbation. As the coupling strength of the network was increased, the PCF reached a maximum at an intermediate coupling strength, while the order parameter increased monotonically (dotted line in Figure 2A). In Figure 2A, we illustrated the intermediate coupling strength (blue region) as the critical state of this brain network model and selected an incoherent (green region) and highly synchronized (red region) states of lower PCFs as a supercritical and subcritical state, respectively. The Φ value was also maximized at a point between sub-and supercritical states, at a similar coupling strength that corresponds to the maximum PCF. To test the robustness of our method, we applied another method for estimating Φ, which is based on Queyranne's submodular optimization algorithm in Kitazono's Φ-toolbox [49]. We found the result was consistent with the maximum Φ at the critical state (See for details). PCFs (in the sub-and supercritical states) are relatively homogeneous and not similar with the structural brain network.  Figure 2A), a salient structured functional connectivity, which resembles the structural brain network, emerges in the brain network. Color bar indicates the node degree.

A Network Mechanism of the Maximal in a Critical State
In a critical state, a network synchronization is balanced by incoherent and synchronous connections through partial phase locking. At a coarse-grained level, the distribution of incoherent and synchronous connections is shaped by the network topology; oscillations at the nodes with dense connections become slower and more synchronous, while oscillations at the nodes with sparse connections are faster and incoherent [16,26]. As a consequence, a coarse-grained functional network resembles a structural network in a critical state [14,55,56]. Figure 3A presents the Spearman correlation coefficients between the node degrees (in the structural brain network) and the PLIs of the 78 nodes (in the functional brain network) as coupling strength increases. The Spearman correlation coefficient is at a maximum in the critical state (blue region), while both the super-and subcritical states (green and red regions, respectively) are associated with smaller correlations. The results imply that the constraint of the structural network on a functional network is maximized in the critical state, with higher degree nodes having a larger PLI. Figure 3B presents the scatter plots of the PLIs versus the node degrees of the 78 nodes in the three states. A large positive correlation coefficient appears only in the critical state ( Figure 3B blue, r = 0.57, p-value < 0.001). However, the large correlation does not mean that the functional network is static. Figure 3C presents the temporal evolution of the Spearman correlation coefficients between the instantaneous phases of alpha oscillations in the functional brain network and the node degrees  Figure 2A), a salient structured functional connectivity, which resembles the structural brain network, emerges in the brain network. Color bar indicates the node degree.
Here, we assumed that the maximum Φ might arise due to the balance between the large variation of functional network configurations and the large constraints of structural network configurations, as a characteristic of a critical state. The disrupted balance in an incoherent or highly synchronized brain network results in a small Φ. Figure 2B presents the functional brain networks based on the phase lag index (PLI). The functional connectivity at the maximum PCF (in the critical state) resembles the structural brain network of DTI, whereas the functional connectivity at the lower PCFs (in the sub-and supercritical states) are relatively homogeneous and not similar with the structural brain network.

A Network Mechanism of the Maximal Φ in a Critical State
In a critical state, a network synchronization is balanced by incoherent and synchronous connections through partial phase locking. At a coarse-grained level, the distribution of incoherent and synchronous connections is shaped by the network topology; oscillations at the nodes with dense connections become slower and more synchronous, while oscillations at the nodes with sparse connections are faster and incoherent [16,26]. As a consequence, a coarse-grained functional network resembles a structural network in a critical state [14,55,56]. temporal change of the phase lead-lag relationships among 78 brain regions upon the structural brain network. The large variation of functional networks is one of the characteristics of a critical state and measured by a large PCF. As a result, the large variation of the functional brain network at a smalltime scale (~seconds; Figure 3C) under the large constraint from the structural brain network ( Figure  3B) in a large time scale (~minutes) in the critical state may be the network condition for the maximal Φ in the brain network.

EEG Study: Correlation Between Criticality, , and Human Consciousness
To investigate the relationship between criticality, integrated information, and level of human consciousness, we compared the PCF and Φ of high-density EEG and the response rate during general anesthesia. The behavioral response rate, which is inversely proportional to the drug concentration, was used as a surrogate for the level of consciousness. In Figure 4A, the conscious resting state and the conscious recovery state have a higher PCF and Φ than in the unconscious states. For the continuous EEG data, the change of PCF correlates with the change of Φ . Moreover, both the measures reflect the response rate during the significant state transition. These results are consistent  Figure 3A presents the Spearman correlation coefficients between the node degrees (in the structural brain network) and the PLIs of the 78 nodes (in the functional brain network) as coupling strength increases. The Spearman correlation coefficient is at a maximum in the critical state (blue region), while both the super-and subcritical states (green and red regions, respectively) are associated with smaller correlations. The results imply that the constraint of the structural network on a functional network is maximized in the critical state, with higher degree nodes having a larger PLI. Figure 3B presents the scatter plots of the PLIs versus the node degrees of the 78 nodes in the three states. A large positive correlation coefficient appears only in the critical state ( Figure 3B blue, r = 0.57, p-value < 0.001). However, the large correlation does not mean that the functional network is static. Figure 3C presents the temporal evolution of the Spearman correlation coefficients between the instantaneous phases of alpha oscillations in the functional brain network and the node degrees in the structural brain network for the three states. The varying correlation coefficients indicate the temporal change of the phase lead-lag relationships among 78 brain regions upon the structural brain network. The large variation of functional networks is one of the characteristics of a critical state and measured by a large PCF. As a result, the large variation of the functional brain network at a small-time scale (~seconds; Figure 3C under the large constraint from the structural brain network ( Figure 3B) in a large time scale (~minutes) in the critical state may be the network condition for the maximal Φ in the brain network.

EEG Study: Correlation Between Criticality, Φ, and Human Consciousness
To investigate the relationship between criticality, integrated information, and level of human consciousness, we compared the PCF and Φ of high-density EEG and the response rate during general anesthesia. The behavioral response rate, which is inversely proportional to the drug concentration, was used as a surrogate for the level of consciousness. In Figure 4A, the conscious resting state and the conscious recovery state have a higher PCF and Φ than in the unconscious states. For the continuous EEG data, the change of PCF correlates with the change of Φ. Moreover, both the measures reflect the response rate during the significant state transition. These results are consistent with the model prediction and empirically demonstrate for the first time a direct relationship between PCF, Φ, and the level of human consciousness. with the model prediction and empirically demonstrate for the first time a direct relationship between PCF, Φ , and the level of human consciousness. The model study also predicted a large correlation of the structural brain network with the functional brain network in consciousness, which corresponds to the functional brain network in a critical state. However, since we recorded only the scalp EEG, we were not able to test the correlation between functional (PLI) and structural (node degree) brain networks. Instead, we compared the PLI networks of the EEG and a simulated EEG. For the simulation of EEG in the conscious state, we first simulated the source signals in the structural brain network in a critical state (the same simulation in Figures 2 and 3) and then projected the source signals onto the scalp (see the Method in details). If the model prediction is correct, both the PLI networks should be largely correlated.
Indeed, we found high Spearman correlation coefficients between the PLI network of the EEG in the conscious states and the PLI network of the simulated EEG in a critical state. This correlation decreases along with the reduced response rate in the unconscious states ( Figure 4B). The high correlation coefficients in conscious states indicate a strong constraint of the structural brain network on the EEG. The Spearman correlation coefficients between the PLI networks of the EEG and the simulated EEG based on the anatomical brain network and critical state. The conscious states (baseline, induction, and recovery) show larger correlations, which imply a stronger constraint of the structural brain network on the EEG in conscious states. As a result, the balance between the large constraint of the structural brain network and the large repertoire (i.e., large PCF) of the functional brain network may be the network basis of the large Φ of EEG in the conscious brains. The PCF and Φ of 64-channel EEG correlate with the response rate (grey area), which was modulated with increasing anesthetic concentrations. (B) The Spearman correlation coefficients between the PLI networks of the EEG and the simulated EEG based on the anatomical brain network and critical state. The conscious states (baseline, induction, and recovery) show larger correlations, which imply a stronger constraint of the structural brain network on the EEG in conscious states. As a result, the balance between the large constraint of the structural brain network and the large repertoire (i.e., large PCF) of the functional brain network may be the network basis of the large Φ of EEG in the conscious brains.
The model study also predicted a large correlation of the structural brain network with the functional brain network in consciousness, which corresponds to the functional brain network in a critical state. However, since we recorded only the scalp EEG, we were not able to test the correlation between functional (PLI) and structural (node degree) brain networks. Instead, we compared the PLI networks of the EEG and a simulated EEG. For the simulation of EEG in the conscious state, we first simulated the source signals in the structural brain network in a critical state (the same simulation in Figures 2 and 3) and then projected the source signals onto the scalp (see the Method in details). If the model prediction is correct, both the PLI networks should be largely correlated.
Indeed, we found high Spearman correlation coefficients between the PLI network of the EEG in the conscious states and the PLI network of the simulated EEG in a critical state. This correlation decreases along with the reduced response rate in the unconscious states ( Figure 4B). The high correlation coefficients in conscious states indicate a strong constraint of the structural brain network on the EEG.

Discussion
Summary of the findings: Previous studies investigated the relationships between criticality and consciousness [12][13][14][15][16][17]19,23,33,57], or integrated information and consciousness [1][2][3][4][5][6]50,58], or criticality and integrated information [19,[59][60][61], respectively. No one has yet studied the relationships between criticality, integrated information, and consciousness at the same time. In this paper, we demonstrated empirically explicit relationships between criticality, integrated information, and consciousness at various levels of human consciousness modulated with a general anesthetic, and also performed a computational model study to understand how the brain at a critical state facilitates the development of a large amount of integrated information at a network level. In the modeling study, we simulated the brain network activities near and far from a critical state and quantified the criticality and integrated information with PCF and Φ. We found that when the brain network reaches a maximal PCF, the Φ is also maximized under the largest constraint of the structural brain network (that is, with the largest correlation between the functional and structural brain networks). To verify the modeling results, we gradually modulated the level of consciousness in humans with a general anesthetic and compared the PCF and Φ of high-density EEG with the behavioral response rate. We found that, in the conscious resting states, the subjects have the highest PCF and Φs when the behavioral response rates are highest. We also showed that the functional brain network of EEGs in conscious states largely correlates with the functional brain network of the simulated EEGs, which were modeled to reflect the underlying structural brain network in a critical state. From both the modeling study and EEG analysis, we proposed that the balance between the large variation of functional brain network (as measured as the large PCF) and the large constraint from the structural brain network in a critical state is a necessary condition to generate the high Φ of conscious states.
Criticality and maximal Φ: In the application to our brain network model, since we did not give any causal relationships among nodes initially, a non-zero Φ may emerge through only the interaction of 78 nodes in the brain network. This raises the following questions. How does a non-zero Φ emerge spontaneously in the brain network? What is the network basis that facilitates the emergence of a non-zero Φ? In this study, we proposed that the phase lead-lag relationship that is shaped by the structural brain network might create conditions that facilitate the emergence of a non-zero Φ in the brain network.
Previous studies have examined how brain network topology modulates the frequencies and phases of local node dynamics, subsequently shaping a pattern of global information flow and functional connectivity. A mathematical relationship between the node degree and the phase of node dynamics was identified analytically [16,26], which enables us to estimate analytically the phase of a node (in a functional network) with only its node degree and local connectivity (in a structural network). This was tested with diverse brain networks (human, monkey, and mouse), which demonstrated that when considering long-term and spatially coarse-grained brain activities (>minutes) the global phase lead-lag relationship and information flow pattern (measured by transfer entropy and Granger causality) were predictable based only on the structural brain network topologies of the three species. In particular, the difference between the "hub node" and "peripheral node" becomes most salient at a critical state. In other words, the structural complexity of overall coupled node dynamics is maximized at a critical state. The differences between frequency and phase among coupled node dynamics naturally give rise to information flow in the brain network, consequently dividing the nodes in the brain network into 'senders' and 'receivers' of information flow. However, in sub-and supercritical states, there is no information flow between node dynamics because of a highly synchronized state (i.e., no difference that would allow for information flow) and incoherent state (i.e., no interaction that would allow for information flow). Therefore, a non-zero Φ may emerge between these two extreme states and be maximized at a balanced state. Extrapolating from results based on the phase lead-lag relationship, emerging a non-zero Φ without external stimuli is likely not random but rather is organized by a mathematical relationship between network structure and dynamics in a critical state [16,26,62,63].
A large variation under a large constraint: In our modeling study and EEG analysis, we found that a large variation of the functional brain network occurs in spite of a strong constraint from the structural brain network. Notably, the large variation and the strong constraint were observed at different time scales. When we investigated the correlation coefficient between the instantaneous phases and the node degrees of the 78 nodes, the correlation coefficient at each time point varied widely at a short time scale (~seconds; Figure 3C). By contrast, the averaged instantaneous phases and the node degrees at a longer time scale (~minutes) have a large positive correlation coefficient (r = 0.57, Figure 3B), which implied a bias toward positive values in the short time scale variation. Our previous EEG study demonstrated that the correlation coefficient between the EEG network and the structural brain network was pronounced when calculated with large temporal windows (>60 seconds) but diminished when calculated with small windows (<5 s). Interestingly, the scale-dependency appeared only in the baseline conscious states, not in the altered states of consciousness such as anesthetized state, minimally conscious state, and unresponsive wakefulness syndrome [23]. The conscious brain can be characterized by diverse repertoires of global functional connectivity on a short time scale. However, when combining all repertoires of small windows in a large window, the different functional connectivity patterns across small windows may be averaged out and only the common functional connectivity patterns that are constrained by the structural brain network remain. Such an averaged functional connectivity in a large window can be simulated by the mean-field method upon a structural brain network. Nevertheless, how the large variance and the large constraint contribute to the large integrated information in conscious states remains elusive. Further studies are required to understand the mechanistic association between the functional variance, structural constraint, and Φ in a network.
Deviations from criticality and Φ: When brain dynamics deviate from a critical state, the capacity to shape global node dynamics into a structure that resembles the network structure is lost. In a supercritical state, there is no functional interaction between nodes. In a subcritical state, node dynamics cannot be structured due to strong functional interactions between nodes that eliminate functional heterogeneity. As a result, the Φ of a brain not in a critical state-regardless of whether the deviation is toward subcriticality or supercriticality-will decrease in proportion to the distance from the critical state.
Many empirical studies support the criticality hypothesis in consciousness by comparing the dynamics of various conscious states with unconscious states (such as sleep, anesthesia, seizure, and coma). They have commonly shown that unconsciousness is associated with a deviation from criticality, which is quantified with various methods (power law, susceptibility, pair correlation function, and correlation between the functional and structural brain networks) [15,22,40,57,[64][65][66][67][68][69][70][71][72]. However, no one has yet associated criticality, human consciousness, and Φ because of the explosive computational demands required for Φ calculations. In a previous study, we introduced a method that statistically estimates the increase/decrease of Φ for a continuous EEG, termed Φ, and showed the applicability to characterizing levels of consciousness with high-density EEG and also enabled us to study a relationship with network criticality [50]. Our current results show that the empirical and computational model studies support the hypothesis that the characteristic network properties in a critical state naturally give rise to a structured, asymmetric phase lead-lag relationship among oscillators and, in turn, maximize Φ in the conscious brain.
Limitations: This study has several limitations. First, Φ does not measure the precise Φ. Instead, we estimated the relative change of Φ, focusing on the increase/decrease along with the change in criticality in the network model and based on empirical EEG during anesthetic-induced unconsciousness. Even so, this relative measure is sufficient for the purpose of this study, which is to examine the relationship between Φ, criticality, and levels of consciousness. Considering the significant differences among the six versions of Φ in the recent comparison [73], we tested another Φ that was introduced recently [49] and found the results were consistent. However, testing other versions of Φ may enhance the reliability of our study. Second, in this study, we derived a relationship between asymmetric phase lead-lag relationships and Φ, but this relationship did not explain how the phase lead-lag relationship generated a causal relationship defined by information partition for Φ. It may require an analytic study to identify conditions to link phase lead-lag and Φ. Third, to find an association between Φ, criticality, and consciousness, we used the subjects' response rate during exposure to a general anesthetic as a surrogate for the level of consciousness. However, it is well known that unresponsiveness does not necessarily correlate with lack of consciousness. Although previously studied in the context of vegetative and minimally conscious states, our research team has recently identified covert consciousness during propofol sedation in which overt motoric responses were absent but brain network responses suggestive of volition were present [74]. Future studies might incorporate EEG data, behavioral responsiveness, and neuroimaging protocols to determine covert consciousness in order to more precisely identify the relationship between critical dynamics, consciousness, and Φ. Fourth, the number of subjects (7) in the empirical data analysis was low. To improve reliability, more subjects may be required in future studies. Testing the relationship level of consciousness, criticality, and Φ with other states of consciousness such as sleep and vegetative states may also be required. Finally, the structural brain network used in the modeling study includes only the cortex. Including subcortical networks such as thalamocortical and hippocampocortical connections, etc. could improve the modeling performance for complex state transitions during general anesthesia.

Conclusions
We demonstrated for the first time an explicit relationship between criticality, integrated information, and human consciousness with computational modeling and EEG analysis. We propose that network criticality, that is, a balanced state between large variation of functional network configurations and strong constraint of structural network configurations, is a network condition for integrated information. Understanding this relationship may open a new way to study diverse states of consciousness situated near to and far from a critical state in terms of integrated information. It may also provide a theoretical foundation for controlling the level of consciousness and integrated information by modulating criticality at a network level.
Supplementary Materials: The following are available online at http://www.mdpi.com/1099-4300/21/10/981/s1, Table S1. Statistical tests of Pair Correlation Function(PCF) and integrated information(Φ) among supercritical, critical, and subcritical states in model, Table S2. Statistical tests of Pair Correlation Function(PCF) and integrated information(Φ) among baseline, induction, anesthesia, and recovery states in experiment. The baseline, induction, anesthesia, and recovery states are the average of 10-min with 0%, 0.4%, 0.8%, 0% drug concentrations, respectively, Figure S1. In order to test the robustness for the relationship between integrated information and criticality (measured by pair correlation function), we tested another integrated information measure, Φ_SI, recently introduced with its toolbox [50]. We compared the two integrated information measures, (A) Φ in Figure 2 and (B) Φ_SI based on Queyrannes's submodular optimization algorithm [48,50]. For Φ_SI, we iterated the simulation 10 times. The error bar indicates the standard deviation at each coupling strength. Both the integrated information measures based on different algorithms demonstrate the maximal Φ values near the critical state. Notably, we did not remove trivial partitions (e.g., partitioning the 78 nodes into 1 node and 77 nodes) in the calculation of Φ_SI, which frequenctly observed with the model data. However, for the prescise calculation of Φ_SI, it should be handled carefully in the future application.
Author Contributions: Hyoungkyu Kim and UnCheol Lee designed the research and wrote the manuscript; Hyoungkyu Kim performed the simulation and data analysis; UnCheol Lee supervised research.
Funding: This research was funded by the National Institutes Health, grant No. R01 GM09857.