Measuring the dynamic balance of integration and segregation underlying consciousness, anesthesia, and sleep in humans

Dataset-1Dataset-1 has been reported in prior publications26,74,75, wherein it was used to test different hypotheses with methodologies distinct from those employed in the current research. The experimental protocol received approval from the Institutional Review Board (IRB) at the University of Michigan, ensuring compliance with all relevant guidelines and regulations. This study involved twenty-six healthy individuals (13 males and 13 females; age range: 19-34 years; all right-handed). Before participating, everyone provided informed consent and was compensated upon completion of their participation in the experiment. To protect participants’ privacy, stringent confidentiality measures were implemented. Each participant was allocated a unique code number at the initial point of contact, which was consistently used as the sole identifier for their behavioral, physiological, and MRI data, ensuring anonymity throughout the research process.The criteria for including and excluding participants were defined as follows: Eligible participants were healthy right-handed individuals classified as American Society of Anesthesiologists physical status 1, aged between 18 and 40, with a body mass index (BMI) of less than 30. Exclusion criteria encompassed any MRI contraindications such as severe obesity, potential pregnancy, claustrophobia, the presence of metal in the body, anxiety disorders, or any cardiopulmonary conditions; a history of neurological, cardiovascular, or respiratory diseases; significant head trauma involving loss of consciousness; any learning or developmental disabilities; a history of sleep apnea or pronounced snoring; any sensory or motor impairments that could affect participation; inability or unwillingness to refrain from alcohol for 24 h before the MRI appointment; gastroesophageal reflux disease; any history of illicit drug use or a positive result on a drug screening test; tattoos located on the head or neck area; known allergies to eggs; detectable intracranial abnormalities observed in T1-weighted MRI images; or any physical discomfort experienced during the fMRI scans.Prior to the study, all participants were required to adhere to an 8-h fasting period. On the day of the experiment, a preoperative evaluation, including a physical examination, was conducted by an attending anesthesiologist. The experiment was continuously overseen by two experienced anesthesiologists who were on-site to monitor vital signs such as spontaneous breathing, end-tidal CO2, heart rate, pulse oximetry, and the electrocardiogram. Arterial pressure was recorded with an MR-compatible monitoring device. Participants received a local anesthetic through a subcutaneous injection of 0.5 mL of 1% lidocaine, followed by the insertion of an intravenous cannula. Supplemental oxygen was also provided at a flow rate of 2 L min−1 via a nasal cannula.We used propofol as the sedative-hypnotic, widely used in fMRI research to investigate the impact of anesthesia on the human brain. Propofol is an attractive drug for this purpose due to its preserved coupling of cerebral metabolism and cerebral blood flow as well as its capacity for precise titration. The drug is thought to work primarily by potentiation of GABA-A receptor-mediated inhibition across the brain. The administration of propofol was meticulously managed through a target-controlled intravenous bolus and a steady infusion regimen. The specific bolus dose, rate of infusion, and duration were established in advance, following the Marsh pharmacokinetic model as integrated into the STANPUMP software (http://opentci.org/code/stanpump). The dosing strategy, combining bolus and continuous infusion, was adjusted at 5-min intervals to achieve the desired concentration at the effect site. The dosage was increased in 0.4 μg mL−1 increments until the subjects displayed no behavioral responses26. The target concentration was sustained for 21.6 ± 10.2 min, after which the infusion ceased, permitting the subjects to recover consciousness naturally.The assessment of behavioral responsiveness involved participants squeezing a rubber ball, a method used to determine the onset of deep sedation induced by propofol, marked by the lack of behavioral response. The BIOPAC MP160 system, coupled with AcqKnowledge software version 5.0, was employed to quantify these responses, measuring the exerted air pressure in mmHg. During the scan, a series of 60 motor response trials were executed, spaced at approximately 90-s intervals. Each trial commenced with the command “action,” prompting participants to squeeze the ball once. Between these motor assessments, participants engaged in various mental imagery exercises, such as envisioning playing tennis, navigating space, or imagining the act of squeezing. Additional details regarding the experimental setup are provided in previous publications26,75.Data collection was performed at Michigan Medicine, University of Michigan, utilizing a 3 T Philips MRI scanner equipped with a standard 32-channel transmit/receive head coil. High-resolution anatomical images were acquired with T1-weighted spoiled gradient recalled echo (SPGR) imaging, employing parameters of 170 sagittal slices, 1.0 mm slice thickness without interslice gaps, a repetition time (TR) of 8.1 s, an echo time (TE) of 3.7 ms, a flip angle of 8°, a field of view (FOV) of 240 mm, and an image matrix of 256 × 256.Functional brain images were obtained using a gradient-echo echo-planar imaging (EPI) sequence, with settings of 28 slices, TR/TE = 800/25 ms by multiband acquisition (MB factor = 4), slice thickness = 4 mm, in-plane resolution = 3.4 × 3.4 mm, FOV = 220 mm, flip angle = 76°, and a 64 × 64 image matrix. Before the upgrade of MRI hardware, six subjects underwent scans with slightly different parameters: 21 slices, TR/TE = 800/25 ms (MB factor = 3), and slice thickness = 6 mm. Inside the scanner, subjects were instructed to stay awake and keep their eyes closed, receiving verbal cues through headphones. The protocol included four fMRI sessions: a 15-min baseline, a 30-min session during and after propofol infusion, and a final 15-min recovery period. In the group analysis (e.g., Fig. 2c), the first and last 5 min of the unresponsive period were excluded from state-wise averaging to rule out any transient effect.Dataset-2The experimental protocol received approval from the IRB at the University of Michigan, ensuring that all procedures complied with relevant guidelines and regulations. The study enrolled thirty healthy individuals (10 males and 20 females; age range: 18-38 years; all right-handed). Each participant provided informed consent before participating and received compensation upon completion of the study. Three participants were excluded due to technical issues with the MRI scanner or excessive movement during scanning. Privacy details and inclusion/exclusion criteria were similar to that of Dataset-1.During the experiment, participants were asked to rest, engage in a behavioral task, or listen to music as outlined below. The anesthetic administration and monitoring were consistent with that of Dataset-1, with several exceptions. The propofol infusion was manually adjusted to progressively reach target effect-site concentrations of 1.5, 2.0, 2.5, and 3.0 μg mL−1. Each target concentration was sustained for a duration of 4 min, allowing us to gradually adjust the anesthetic dosage to pinpoint the exact threshold for LOR. To reduce the potential for head motion-related artifacts, we kept the ESC at one level higher than the concentration associated with LOR for ≈ 32 min. For instance, if a participant exhibited LOR at a concentration of 2.0 μg mL−1, the effect-site concentration was maintained at 2.5 μg mL−1 for the duration of the LOR phase. In rare cases (occurrence rate of ~5.6%) where participants did not lose their responsiveness at 3.0 μg mL−1, the target concentration was escalated to 3.5 μg mL−1 and sustained at maximum of 4.0 μg mL−1. Following this, the infusion was discontinued, allowing the propofol levels to diminish gradually.Throughout the experiment, a series of eight fMRI scans were performed, each lasting 16 min. These scans comprised a 16-min resting-state (Rest1) and a 16-min music-listening (Music1) session during the baseline period of wakefulness, followed by a 16-min behavioral test conducted during the induction phase of propofol. After the participants had lost behavioral responsiveness, they underwent another set of 16-min resting-state (Rest2) and music-listening (Music2) scans. This was succeeded by a 16-min behavioral test during the emergence phase, which occurred after termination of the propofol infusion. The final scans included a 16-min resting-state (Rest3) and a 16-min music-listening (Music3) session after the participants had recovered behavioral responsiveness. Breaks ranging from 1 to 5 min were present between each scan. The entire imaging protocol and data collection process were completed within a total scanning session of 2.5 h for each participant. In the group analysis (e.g., Fig. 2d), Rest1 and Music1 were defined as baseline state. Rest2 and Music2 were defined as LOR state. Rest3 and Music3 were defined as recovery state.In the resting-state period, participants were asked to lie still with their eyes closed, trying to remain awake while refraining from engaging in specific thoughts or movements. For the music periods, they were asked to listen attentively to well-known tracks from four music genres—Jazz, Rock, Pop, and Country—with their eyes closed, remaining motionless and awake. The music tracks were edited to a consistent length of 4 min and played in a pseudo-random sequence. During the behavioral testing periods, participants were prompted to squeeze an MRI-compatible grip dynamometer (a rubber ball) at 10-s intervals for a total of 96 cycles. The cue to initiate each squeezing cycle was the spoken command “squeeze,” delivered through an MRI-compatible audiovisual system, with programming by E-Prime 3.0 software (Psychology Software Tools, Pittsburgh, PA). Headphone volume was set to ensure comfort during wakefulness and increased to 150% following the loss of responsiveness to ensure audibility. The participants’ grip strength was quantified in mmHg using the BIOPAC MP160 system coupled with AcqKnowledge software (V5.0). Analysis of the synchronization between the “squeeze” cues and the participants’ actual hand movements throughout the graded propofol administration allowed for the identification of the moments when subjects transitioned into and out of unresponsiveness. These transitions were marked by the first missed and successfully completed hand squeeze, respectively.Data collection was conducted at Michigan Medicine, University of Michigan, using a 3 T Philips MRI scanner equipped with a standard 32-channel transmit/receive head coil. High-resolution anatomical images were captured using T1-weighted SPGR imaging techniques, configured with the following parameters: 170 sagittal slices, a slice thickness of 1.0 mm, a TR of 8.1 s, a TE of 3.7 ms, a flip angle of 8°, a FOV of 240 mm, and an image matrix size of 256 × 256. Functional brain imaging was performed employing a gradient-echo EPI sequence, detailed by the following specifications: 40 slices, TR/TE = 1400/30 ms by multiband acquisition (MB factor = 4), a slice thickness of 2.9 mm, an in-plane resolution of 2.75 × 2.75 mm, FOV of 220 mm, a flip angle of 76°, and an 80 × 80 image matrix.Dataset-3 (light and deep sedation)This dataset, of which the protocol approved by the IRB of Medical College of Wisconsin (MCW), has been previously published with analyses that differ from those employed in this study76,77. The study involved fifteen healthy participants (9 males and 6 females; age range: 19-35 years; all right-handed) who underwent sedation with propofol. The Observer’s Assessment of Alertness/Sedation (OAAS) scale was utilized to assess behavioral responsiveness. Participants were fully responsive to verbal stimuli with an OAAS score of 5 during the baseline and recovery conditions. In the light sedation condition, they exhibited lethargic responses to verbal commands, marked by an OAAS score of 4, and during deep sedation, they did not respond, with OAAS score ranging from 1 to 2. The target plasma concentrations for propofol varied among individuals (light sedation: 0.98 ± 0.18 μg mL−1; deep sedation: 1.88 ± 0.24 μg mL−1), reflecting the inter-individual differences in sensitivity to the anesthetic. To maintain a steady state of sedation, the propofol plasma concentration was regulated by adjusting the infusion rate, balancing the drug’s accumulation and elimination, with manual adjustments informed by STANPUMP. Standard American Society of Anesthesiologists monitoring was conducted, which included tracking the electrocardiogram, blood pressure, pulse oximetry, and end-tidal CO2, along with prophylactic supplemental oxygen delivered through a nasal cannula. The study had to exclude one participant due to excessive movement, leaving 14 subjects available for subsequent analysis.The acquisition of resting-state data encompassed four 15-min scans, each corresponding to different conditions: baseline consciousness, light sedation, deep sedation, and recovery. The images were captured using a 3 T Signa GE 750 scanner (GE Healthcare, Waukesha, Wisconsin, USA) equipped with a standard 32-channel transmit/receive head coil. The scanner obtained whole-brain gradient-echo EPI sequences, detailed as follows: 41 slices per brain volume, TR = 2 s, TE = 25 ms, slice thickness = 3.5 mm, FOV = 224 mm, flip angle = 77°, and an image matrix size of 64 × 64. Alongside these functional scans, high-resolution anatomical images were captured for coregistration.Dataset-4 (light sedation and surgical level of propofol)The dataset has been published previously, utilizing analyses that differ from the ones applied in the current study78. The study was approved by the IRB of Huashan Hospital, Fudan University. The study included twenty-six individuals (12 males and 14 females; age range: 27-64 years; all right-handed). Prior to their participation, informed consent was secured from all candidates, who subsequently received compensation post-experiment. All participants were classified as American Society of Anesthesiologists (ASA) physical status I or II and were scheduled for elective surgery via a trans-sphenoidal approach for the removal of pituitary microadenomas. The diagnosis of these microadenomas was confirmed through radiological evaluations and plasma endocrine levels, identifying tumors less than 10 mm in diameter that did not extend beyond the sellar region. Individuals with history of brain dysfunctions, significant organ failure, usage of neuropsychiatric medications, MRI contraindications like metal implants or vascular clips were excluded. Three individuals were subsequently excluded from the analysis due to excessive movement, leaving data from 23 participants for further examination.Prior to the study, participants were instructed to abstain from solid foods for at least 8 h and from liquids for a minimum of 2 h. Throughout the fMRI session, continuous monitoring was conducted for blood pressure, electrocardiography, pulse oximetry, and the partial pressure of CO2. An intravenous catheter was inserted into a vein on the right hand or forearm, through which propofol was administered to 17 out of 23 participants for light sedation and to all 23 participants for general anesthesia. Propofol administration was managed through a target-controlled infusion system, calibrated according to the Marsh model to sustain a stable ESC. For general anesthesia procedures, remifentanil was administered at 1.0 μg kg−1, along with succinylcholine at 1.5 mg kg−1, to aid in endotracheal intubation. The infusion began at a starting concentration of 1.0 μg mL−1, with 0.1 μg mL−1 increments applied until the target ESC was reached. A 5-min interval was allowed for the propofol to evenly distribute across the compartments. The target propofol ESCs were maintained at 1.3 μg mL−1 for light sedation and 4.0 μg mL−1 for surgical levels of general anesthesia.Ramsay scale was used to evaluate the behavioral responsiveness of the participants. Participants were categorized as fully conscious (Ramsay scores 1-2) when they provided a clear and prompt response to the command (“strongly squeeze my hand!”). If their response was clear but slow, they were considered to be mildly sedated (Ramsay scores 3-4). Those who did not respond were classified as deeply sedated or anesthetized (Ramsay scores 5–6). This specific verbal command of the Ramsay scale was administered twice during each evaluation phase. In conditions of conscious resting state and light sedation, participants were allowed to breathe spontaneously, receiving supplemental oxygen through a nasal cannula. During general anesthesia, ventilation was supported with intermittent positive pressure, maintaining a tidal volume of 8-10 mL kg−1 and a respiratory rate of 10-12 breaths per minute. The procedure was overseen by two certified anesthesiologists throughout its duration. To ensure auditory isolation and comfort during the fMRI scans, participants were equipped with earplugs and headphones.Three fMRI scans (each 8-min long) were performed across different states: conscious baseline, light sedation, and general anesthesia. To minimize the head movement, participants’ heads were securely positioned within the scanner frame and supported by soft padding. They were instructed to lie back comfortably, remain supine, and keep their eyes closed, with an eye patch applied. During the resting-state scans, participants were advised to relax and not focus on any specific thoughts. The imaging was performed using a Siemens 3 T MAGNETOM scanner equipped with a standard 8-channel head coil, capturing whole-brain gradient-echo EPI images under the following settings: 33 slices, a TR of 2000 ms, a TE of 30 ms, a 5 mm slice thickness, a 210 mm FOV, a 90° flip angle, and a 64 × 64 image matrix. High-resolution anatomical images were also obtained. For the current study, we used the fMRI data collected during the conscious baseline and general anesthesia.Dataset-5: Natural sleepThe dataset was acquired from the OpenNEURO database, shared by Pennsylvania State University with informed consent40,41. The original dataset included 33 healthy participants with wakefulness and three sleep stages (N1, N2, and N3). The experimental procedures were approved by the IRB at Pennsylvania State University. Sleep stage was identified based on EEG signatures determined by a Registered Polysomnographic Technologist40,41. Anatomical scans were conducted using a 3 Tesla Prisma Siemens Fit scanner, equipped with a Siemens 20-channel receive-array coil. Functional scans were performed through an EPI sequence, with a TR of 2.1 s, a TE of 25 ms, a slice thickness of 4 mm, comprising 35 slices, a FOV of 240 mm, and an in-plane resolution of 3 × 3 mm. One subject with head motion greater than 10 mm was excluded from the analysis. Only three out of 33 participants exhibited N3 sleep stages in their fMRI data, making analysis of this stage unreliable and leading to its exclusion. For ISD analysis, we excluded an additional four subjects who had any of the stages (awake, N1, and N2) shorter than 1 min (n = 28). For metastability and complexity analysis, we excluded 12 subjects who had any of the stages (awake, N1, and N2) shorter than 5 min (n = 20).Dataset-HCPThe dataset was acquired from the S1200 Release of the WU-Minn Human Connectome Project (HCP) database, which has been comprehensively described in previous publications79. Participants were healthy young adults from 22 to 37 years old. Participants who completed two sessions of resting-state fMRI scans (Rest1 and Rest2) were selected (n = 1009). Data collection utilized multiband EPI via a customized Siemens 3 T MR scanner (Skyra system). Each scanning session included two sequences, each with different phase encoding directions (left-to-right and right-to-left), lasting 14 min and 33 s, with a TR of 720 ms, a TE of 33.1 ms, and a voxel size of 2-mm isotropic. The sequences from each session were merged, totaling 29 min and 6 s of data per session. This approach of merging data from opposite phase encodings aimed to neutralize any biases introduced by the direction of phase encoding. The denoised volumetric data, processed through ICA-FIX, were obtained from the online HCP database. Further details on the data collection and preprocessing procedures are documented in other publications80,81. Static functional connectivity matrices of the Rest1 period produced by Pearson correlation were used for the analysis.Sex and gender informationFor Dataset-1 to Dataset-4, participant sex was determined by self-report, acknowledging that self-reported sex may not always reflect biological sex. Gender identity was not collected in this study. While our design included both male and female participants, sex-based analyses were not planned in advance due to the nature of the study design and the limited sample size, potentially precluding meaningful post-hoc analysis. However, Source Data is reported disaggregated by sex to promote transparency and facilitate future investigation. For Dataset-5, individual-level demographic information was not available. For Dataset-HCP, participants self-identified their biological sex as male or female. Gender identity was not assessed. (https://www.humanconnectome.org/storage/app/media/documentation/s1200/HCP_S1200_Release_Reference_Manual.pdf).fMRI data preprocessing and regional time course extractionThe preprocessing of fMRI data adhered to a comprehensive protocol established in our previous studies26,74,82. The AFNI software suite (linux_ubuntu_16_64; http://afni.nimh.nih.gov/) was utilized. Steps included: (1) Slice timing correction; (2) Rigid body head motion correction and realignment; (3) Frame-wise scrubbing to address head motion artifacts; (4) Coregistration with T1-weighted anatomical images; (5) Spatial normalization into Talairach stereotactic space, resampled to 3 × 3 × 3 mm; (6) Band-pass filtering to 0.01-0.1 Hz while regressing out unwanted noise such as drifts or motion artifacts, and physiological noise from white matter and cerebrospinal fluid signals; (7) Spatial smoothing with 6 mm full width at half maximum isotropic Gaussian kernel; (8) Temporal normalization to zero mean and unit variance. Cortex was defined by 400 cortical areas according to a well-established cortical parcellation scheme83,84. Subcortex was defined by 50 areas derived from a Subcortical Atlas85. The fMRI time courses were extracted from those 450 regions of interest (ROIs).Sliding window analysisWe adopted a sliding window approach for dynamic analysis. We chose a window size of 4 min with a 1-TR time step to reflect the dynamic changes in brain activity caused by the anesthetics instead of intrinsic variabilities73. However, we verified that the performance of ISD in discriminating awake and LOR states was stable across a range of window sizes from 0.5 to 8 min (Fig. 3d and Supplementary Fig. 5d). Pearson correlation was used to construct functional connectivity matrices.Global signal regression (GSR)To assess different properties of functional brain networks, we utilized global signal regression (GSR). GSR is an effective way to mitigate confounding effects of motion that can affect the accuracy of network topology measurements. Moreover, since GSR artificially removes the global temporal coordination of fMRI signals, it can act as an effective way to regress the ‘integration’ of the network and reveal the underlying topological architecture (e.g., segregation), independent of global connectivity. Therefore, we utilized BOLD data that had not been subjected to GSR (noGSR) for calculating integration. For the computation of segregation, we utilized GSR-applied data.Multi-level efficiency and multi-level clustering coefficientWe propose multi-level efficiency and multi-level clustering coefficient as proxies for network integration and segregation, respectively. Refer to Supplementary Fig. 12 for an overview of the calculation procedures. A functional connectivity matrix was binarized at a given threshold t ranging from 0 to 1 (with an interval of 0.01). Global efficiency was calculated from the binary matrix:$${E}_{{glob}}(t)=\frac{1}{N(N-1)}{\sum}_{i\ne j}^{N}\frac{1}{{d}_{{ij}}\left(t\right)}$$
(1)
where \({d}_{{ij}}\left(t\right)\) is the shortest path length between nodes i and j when binarized at threshold t, and N is the total number of ROIs (= 450). Also, global clustering coefficient was obtained as a mean of the local clustering coefficients:$${C}_{i}\left(t\right)=\frac{2{e}_{i}\left(t\right)}{{k}_{i}\left(t\right)\cdot ({k}_{i}\left(t\right)-1)}$$
(2)
$${C}_{{glob}}\left(t\right)=\frac{1}{N}\sum {C}_{i}\left(t\right)$$
(3)
where \({k}_{i}\left(t\right)\) and \({e}_{i}\left(t\right)\) are the degree of node i and the number of edges between the neighbors of node i when binarized at threshold t.Multi-level efficiency (i.e., integration) and multi-level clustering coefficient (i.e., segregation) were determined by calculating the areas under the curve of the global efficiency and clustering coefficient values across the threshold range.$$\left({Integration}\right)={\int }_{\!\!\!0}^{1}{E}_{{glob}}\left(t\right){dt}$$
(4)
$$\left({Segregation}\right)={\int }_{\!\!\!0}^{1}{C}_{{glob}}\left(t\right){dt}$$
(5)
This method allows for a nuanced analysis of functional connectivity matrices by assessing network properties at various thresholds, thereby capturing the complexity of brain networks beyond single-scale approaches. Lastly, ISD was defined as the multi-level efficiency subtracted by multi-level clustering coefficient. This multi-thresholding approach is comparable with a previously suggested method of clustering coefficient calculation for complete weighted network44.As both global efficiency and clustering coefficient are affected by the density of functional connectivity matrix, we evaluated possible correlations of integration and segregation with connectivity values (i.e., mean value of positive edges) of matrices with and without GSR (Supplementary Fig. 13). High Spearman correlations were observed between integration and noGSR connectivity, and between segregation and GSR connectivity. This agrees with the intuitive concept of integration and segregation, as integration can be understood as the amount of common global signal generated by a neural system. After GSR, the remaining signal likely reflects more localized, within-module connectivity18, as the global signal has been removed. Increased connectivity after GSR suggests enhanced strength of these local, segregated networks, indicating a shift towards greater segregation. Therefore, this demonstrates that our metrics effectively captured the integration and segregation of a given network.Brain network visualizationBrain plots were generated using BrainNet Viewer (https://www.nitrc.org/projects/bnv), in which connections stronger than r = 0.75 are shown. Topological profile of brain networks was visualized with the open-source software Gephi 0.10.1 (https://gephi.org/), in which the top 10% of connections are depicted using the ForceAtlas2 algorithm.Modularity and participation coefficientModularity is a measure that characterizes how well a network can be divided into modules, which is a widely used metric for segregation. Modularity was defined by:$$Q=\frac{1}{2E}{\sum}_{i,j}\left({A}_{{ij}}-\frac{{k}_{i}{k}_{j}}{2m}\right){\delta }_{{ij}}$$
(6)
where E is total number of edges and \({k}_{i}\) is the degree of node i. The \({\delta }_{{ij}}\) is 1 if i and j belong to the same network, and otherwise 0. For module assignment, eight predefined networks were used.Participation coefficient quantifies between-module strength and has been widely used to measure network integration16,18,28. The participation coefficient of each node was calculated as:$${{PC}}_{i}=1-{\sum}_{m\in M}{\left(\frac{{k}_{i}^{m}}{{k}_{i}}\right)}^{2}$$
(7)
where m is the network in which i does not belong, and \({k}_{i}^{m}\) is the sum of connections between node i and nodes in system m18. The average value of PCi across the nodes was then used as a proxy for network integration.Cartographic profilingCartographic profiling involves constructing 2D joint histograms of the participation coefficient and within-module z-score from every time frame and clustering these histogram profiles into two clusters, i.e., predominantly integrated and segregated states28,86. We followed the same procedure while using predefined 8-network assignment for participation coefficient and within-module z-score calculation. The k-means clustering (k = 2) was performed for each subject and was replicated 500 times with random initialization. The cluster with the higher average participation coefficient was assigned as the predominantly integrated state, while another cluster was labeled as the predominantly segregated state.System segregationSystem segregation quantifies the strength of within-system connections in relation to between-system connections18. It is defined as$$\left({System\; Segregation}\right)=\frac{{\bar{z}}_{{within}}-{\bar{z}}_{{between}}}{{\bar{z}}_{{within}}}$$
(8)
where \({\bar{z}}_{{within}}\) is the mean value of intra-network Fisher z-transformed correlation values and \({\bar{z}}_{{between}}\) is the mean value of inter-network Fisher z-transformed correlation values. Following the original definition18, negative weights were excluded from calculation. Each node was assigned to a system based on the eight predefined networks.Estimation of temporal regimes of ISD changesTo increase the statistical power, Datasets 1 and 2 were combined (total n = 45). For induction, the second scan from Dataset-1 and third scan from Dataset-2 were used. For emergence, the third scan from Dataset-1 and sixth scan from Dataset-2 were used. Each segment was time-locked to the onset of transition. If the lengths of the data before or after transition were shorter than 8 min, the empty values were filled with the mean ISD trends to ensure continuity of smoothing. Each curve was then spline-smoothed with the smoothing parameter of 10-6, which reasonably preserves the general trend while suppressing transient fluctuation. The slope was estimated by obtaining the time derivative of the smoothed curves. For each second, the distribution of 45 slope values was tested against zero median (Wilcoxon signed-rank test, FDR-corrected p < 0.05). Log10 of p-values were applied to moving average of a 4-min window.The goodness of fit was quantified as normalized mean absolute error (NMAE), defined as:$${NMA}{E}_{i}=\frac{1}{n}{\sum}_{t}\frac{\left|{S}_{i}\left[t\right]-\bar{S}\left[t\right]\right|}{\max \left({S}_{i}\right)-\min ({S}_{i})}\times \left(100\%\right)$$
(9)
where \({S}_{i}\left[t\right]\) is spline-smoothed ISD curve of ith subject and \(\bar{S}[t]\) is the group average of spline smoothed curves.Training and evaluation of machine learning modelsIntegration and segregation values were extracted from the dynamic functional connectivity matrices of nine networks (whole-brain and eight pre-defined networks). For model comparison, we aggregated the integration and segregation across all networks, culminating in a feature set of 18 measurements. For training, data points from stable awake and LOR states (Dataset 1: Baseline, Deep LOR [LOR excluding first and last 5 min], and Recovery; Dataset 2: Rest1, Music1, Music3, and Rest3) were used. Data points close to state transitions were not provided to the model during training. Considering the class imbalance in the training set, non-overlapping 20,000 data points were randomly selected from each class (awake and LOR) during training for all five models (total 40,000 training data points).The performance of models was evaluated by the AUC of the ROC curve and balanced accuracy (the mean of sensitivity and specificity). We used a leave-one-subject-out strategy for cross-validation and hyperparameter selection87. This method involved sequentially excluding one subject from the training process and utilizing the omitted subject’s data as the validation set. To ensure the generalizability of the model, hyperparameters chosen from Dataset-1 were also applied to Dataset-2.For random forest, the number of trees is set at 100. Trees are allowed to grow until all leaves are pure. For each split, five features are randomly sampled, equivalent to the ceiling of the square root of the number of total features. For support vector machine, we selected a radial basis function kernel. The regularization parameter C was set to 1.0. Kernel scale parameter was heuristically chosen by fitcsvm function in MATLAB. For artificial neural network, we implemented a pattern recognition network with two hidden layers (10 and 2 neurons each), utilizing a scaled conjugate gradient training function and cross-entropy as the performance function. The learning rate and momentum was set as 0.01 and 0.9, respectively. We employed L2 regularization to prevent overfitting. Our artificial neural network was configured to use rectified linear units for activation and included dropout and early stopping to enhance performance. For k-nearest neighbor model, we selected the optimal value of k = 5 via cross-validation. We employed the Euclidean distance metric due to its effectiveness in continuous data spaces. For logistic regression, we employ a logit link function.MetastabilityWe calculated the moving standard deviation of Kuramoto order parameters as a proxy to quantify metastability in the neural system31. Since GSR removes the entire global synchronization, the noGSR data were used for metastability calculation. The BOLD time series of ROIs were transformed to coupled oscillators via Hilbert transform. The Kuramoto order parameter (R) was defined as:$$R\left(t\right)=\frac{1}{N}\left|{\sum}_{j=1}^{N}{e}^{i{\theta }_{j}\left(t\right)}\right|$$
(10)
where \({\theta }_{j}\) is the instantaneous phase of jth ROI. Metastability was quantified as the moving standard of the Kuramoto order parameter (\({\sigma }_{R}\)) with the sliding window of 4 min.Pattern complexityWe introduce a metric to quantify the complexity of the BOLD signal pattern. The parcellated BOLD signal can be considered as 450-dimensional vectors, which can be grouped by the k-means clustering method, a commonly used unsupervised machine learning algorithm. The clustering was performed at the subject-level, in which the value of k ranged from 3 to 100. To ensure convergence, the clustering was replicated 100 times for each k. Considering the data’s high dimensionality, correlation distance was used as a distance metric. This converts the BOLD signal into a pattern sequence.Subsequently, we calculated three different aspects of complexity with the sliding window of 4 min. First, the entropy of pattern occurrence was calculated as:$${H}_{{patt}}(k)=\frac{-{\sum }_{i=1}^{k}{p}_{i}{\log }_{2}{p}_{i}}{{\log }_{2}k}$$
(11)
where \({p}_{i}\) is the occurrence rate of ith pattern within the given sequence. This quantifies how evenly every pattern appears. Second, we transformed the pattern sequence into k-by-k transition probability matrix and calculated its entropy as:$${H}_{{TM}}(k)=\frac{-{\sum }_{i,j=1}^{k}{p}_{{ij}}{\log }_{2}\,{p}_{{ij}}}{{\log }_{2}\left({k}^{2}\right)}$$
(12)
where \({p}_{{ij}}\) is the occurrence rate of i-to-j transition. Third, we used an effort-to-compress (ETC) technique to quantify the temporal complexity of the pattern sequence37,39. ETC utilizes a lossless compression through non-sequential recursive pair substitution and describes how many substitutions are needed to fully describe the categorical sequence data. Following the consideration of the previous study37, we normalized the raw ETC by dividing with \(L-1\), where L is the length of the series. We combined them by averaging their normalized AUC. We named this composite measure as pattern complexity (Cpatt).$${C}_{{patt}}= \frac{1}{3}\left({nAU}{C}_{{H}_{{patt}}}+{nAU}{C}_{{H}_{{TM}}}+{nAU}{C}_{{ETC}}\right) \\= \frac{1}{3}\cdot \frac{1}{98}{\sum }_{k=3}^{100}\left({H}_{{patt}}\left(k\right)+{H}_{{TM}}\left(k\right)+\frac{{ETC}\left(k\right)}{L-1}\right)$$
(13)
The strong mutual correlations among the entropy of pattern occurrence (Hpatt), transition probability matrix (HTM), and ETC (Supplementary Fig. 11) suggest that our pattern complexity measure (Cpatt) is a robust indicator of the brain’s dynamical state. We confirmed that GSR had a negligible effect on the pattern complexity.Statistics and reproducibilityFor all statistical analyses in Figs. 2d–f, 6f, 7e, 8a, b, d, measurements during three conditions (baseline, LOR, and recovery for anesthesia data; awake, N1, and N2 sleep for sleep data) were statistically compared using Friedman’s ANOVA and post-hoc Wilcoxon test (two-sided). The statistical significance was set at FDR-corrected p < 0.05. In Fig. 4, 45 slope values from each individual for a given time point were tested against zero median by Wilcoxon signed-rank test (two-sided). Multiple comparison was conducted for all networks and all time points (FDR-corrected p < 0.05). In Supplementary Fig. 6, six groups were compared using Kruskal-Wallis test and post-hoc Dunn’s test (two-sided).In Figs. 6i, 7h, 8e,g, we employed a dominance analysis35 to investigate the relationship between metastability/complexity (the response variables) and two key predictors: integration and segregation. Dominance analysis quantifies the relative importance of predictors in a multiple linear regression model by examining their contribution to the model’s explanatory power (adjusted R2). To achieve this, the model is fitted to all possible combinations of predictors, and the average increase in R2 resulting from adding a specific predictor to each submodel is calculated. This average increase represents the predictor’s total dominance. The sum of all predictors’ dominance values equals the total adjusted R2 of the full model, offering a clear way of attributing the model’s overall performance to individual predictors.While individual fMRI data underwent standard motion correction (e.g., motion correction, scrubbing, regressing out motion profiles), we conducted additional group-level control analyses to account for subject- or state-specific motion effects on ISD, metastability, and pattern complexity. Linear regression, using state-averaged FD values as a second-level regressor, isolated fMRI-derived measures from possible motion effects. ANOVA and post-hoc tests on these adjusted values were performed to examine the robustness of our findings (Source Data).The key finding of the paper, i.e., decrease of ISD under LOR induced by propofol or natural sleep, was first found in Dataset-1, which was reproduced independently in Dataset-2 through Dataset-5. Another key finding of the paper, near-zero value of ISD of healthy awake brain was first found in Dataset-1 and replicated independently in Dataset-2, Dataset-5, and Dataset-HCP. Relationships among metastability, complexity, integration, and segregation were first found in Dataset-1 and independently replicated in Dataset-2 and Dataset-5.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.