Muting, not fragmentation, of functional brain networks under general anesthesia

Changes in resting-state functional connectivity (rs-FC) under general anesthesia have been widely studied with the goal of identifying neural signatures of consciousness. This work has commonly revealed an apparent fragmentation of whole-brain network structure during unconsciousness, which has been interpreted as reflecting a break-down in connectivity and disruption in the brains ability to integrate information. Here we show, by studying rs-FC under varying depths of isoflurane-induced anesthesia in nonhuman primates, that this apparent fragmentation, rather than reflecting an actual change in network structure, can be simply explained as the result of a global reduction in FC. Specifically, by comparing the actual FC data to surrogate data sets that we derived to test competing hypotheses of how FC changes as a function of dose, we found that increases in whole-brain modularity and the number of network communities considered hallmarks of fragmentation are artifacts of constructing FC networks by thresholding based on correlation magnitude. Taken together, our findings suggest that deepening levels of unconsciousness are instead associated with the increasingly muted expression of functional networks, an observation that constrains current interpretations as to how anesthesia-induced FC changes map onto existing neurobiological theories of consciousness.


Introduction
While much work has focused on how different anesthetics affect ion channels and receptor function at the cellular level ( Anis et al., 1983;Franks, 2006;Peduto et al., 1991 ), it remains poorly understood, by comparison, how anesthetics affect the coordinated activity of distributed whole-brain networks ( Alkire et al., 2008;Brown et al., 2011 ). In recent years, resting-state functional MRI (rs-fMRI) has provided important glimpses into the large-scale, network-level effects of anesthesia. This approach, which measures covariance structure in spontaneous low-frequency oscillations in neural activity ( Biswal et al., 1995 ), has repeatedly revealed an apparent fragmentation of functional brain network structure during various states of unconsciousness ( Boly et al., 2012b;Hutchison et al., 2014 ). These findings are broadly consistent with work using electroencephalography and electrocorticography reporting a similar breakdown in long-range corti-accounts (e.g. Tononi, 2004 ), the fragmentation of network structure observed during sleep or anesthesia can be viewed as a causal signature of unconsciousness, resulting in a disruption of the brain's ability to integrate information, or to broadcast it widely enough to create conscious awareness ( Mashour, 2013;Mashour and Hudetz, 2018 ).
Work attempting to relate features of functional brain networks to consciousness faces two broad challenges: First, it must faithfully recover and summarize the underlying network structure which gives rise to the observed brain activity; and second, it must determine which components of this structure are causally related to conscious experience, and are not merely epiphenomenal. It is important to distinguish between these two broad goals, given that a misleading characterization of brain network structure necessarily undercuts any interpretations, consciousness related or otherwise, that stem from that network characterization. It is this first goal that the current study seeks to address; namely, deriving a faithful characterization of the whole-brain network changes that underlie the patterns of brain activity commonly observed across neuroimaging studies of consciousness/unconsciousness. Specifically, we investigate the wellcharacterized phenomena of network fragmentation, and the extent to which this is an accurate characterization of changes in network activity under anesthesia.
Network fragmentation is often concluded on the basis of changes in graph properties (e.g., increases in modularity or the number of communities) estimated from thresholded correlation matrices using a fixed, magnitude threshold. However, in this regime, global decreases -or "muting " -of functional connectivity may also produce the appearance of fragmentation by producing sparser, more disconnected networks (see Fig. 1 a). As an illustrative example, we contrast work which observed an increase in the modularity and number of communities of significance thresholded networks during sleep ( Boly et al., 2012b ), with work which found no such effect when using a relative threshold, chosen to obtain a fixed edge density ( Uehara et al., 2014 ). It is diffcult to reconcile these results, or to establish that they represent truly distinct phenomena, as doing so requires a detailed analysis of the underlying network structure. This consideration is particularly pertinent in the case of anesthesia, which has been observed to result in overall decreases in correlation magnitude ( Bettinardi et al., 2015;Lv et al., 2016 ). Further, while several authors have highlighted diminished long-range connectivity under anesthesia -notably between the frontal and parietal cortices  ) -other work has observed similar decreases in local connectivity ( Monti et al., 2013 ). This suggests that anesthetic compounds may affect both short and long range cortical connections, and thus do not simply result in the disconnection of distant cortical regions, but disrupt coordinated neural activity at the local level Vizuete et al., 2014 ).
For studies seeking to identify the neural correlates of unconsciousness using functional connectivity, these facts suggest an important distinction -between alterations in network structure on one hand, versus overall changes in correlation magnitude on the other. These different effects can be difficult to disentangle, as evidenced by findings of both (1) an increase in the modularity and number of communities during sleep, and (2) that, despite these network changes, overall network structure remained relatively preserved ( Boly et al., 2012b ). Given this general ambiguity, as well as the ambiguity as to how such results may support theories of consciousness ( Tononi, 2004 ), the goals of the present study were two-fold: First, to ascertain to what extent network fragmentation under anesthesia-induced unconsciousness is attributable to an overall decrease in connectivity strength; and second, to characterize the structure of whole-brain functional connectivity across depths of unconsciousness. To unpack these relationships, we examined changes in rs-fMRI brain network structure in non-human primates across six increasing levels of anesthesia. This allowed us to test for fine-graded changes in network strength and structure across sedation levels, while also assessing critical components of neurobiological theories related to consciousness.

Data collection
We reanalyzed data from five Macaque primates (M. Fascicularis; 4 female; Mean age 7.8 yrs) collected as part of the experiment reported in Hutchison et al. (2014) . All surgical and experimental procedures were carried out in accordance with the Canadian Council of Animal Care policy on the use of laboratory animals and approved by the Animal Use Subcommittee of the University of Western Ontario Council on Animal Care. As data acquisition is thoroughly described by the original authors, we present a more condensed description here. Note however, that our fMRI preprocessing pipeline contains slight differences.
Prior to image acquisition, subjects were injected intramuscularly with atropine (0.4 mg/kg), ipratropium (0.025 mg/kg), and ketamine hydrochloride (7.5 mg/kg), followed by intravenous administration of 3 ml propofol (10 mg/ml) via the saphenous vein. Subjects were then intubated and switched to 1.5% isoflurane mixed with medical air. Each subject was then placed in a custom-built chair and inserted into the magnet bore, at which time the isoflurane level was lowered to 1.00%. Prior to image localization, shimming, and echo-planar imaging (EPI), at least 30 min was allowed for the isoflurane level and global hemodynamics to stabilize at this concentration. We then acquired 2 functional EPI scans at each of six increasing isoflurane levels: 1.00, 1.25, 1.50, 1.75, 2.00, and 2.75% (0.78, 0.98, 1.17, 1.37, 1.56, and 2.15 minimum alveolar concentration, respectively). We interleaved a 10 min period between each isoflurane level increase to allow for the concentration to stabilize, during which no fMRI data were collected. Throughout the duration of scanning, the monkeys spontaneously ventilated and we monitored physiological parameters (temperature, oxygen saturation, heart rate, respiration, and end-tidal CO2; see Supplemental Figure S2) to ensure that values were within normal limits. The acquisitions of two anatomical images occurred during the stabilization periods between isoflurane levels.

fMRI preprocessing
Functional image preprocessing was implemented using Nipype (Neuroimaging in Python: Pipelines and Interfaces; http://nipy.org/nipype ). The functional images underwent de-spiking, motion correction (six-parameter affine transformation), and slice-time correction, before brain extraction, and highpass temporal filtering (0.01 Hz; Hallquist et al., 2013;Power et al., 2014 ). Functional data was co-registered to its respective T1 anatomical (six degrees of freedom rigid transformation), and then linear (12 degrees of freedom linear affine transformation) and nonlinear transformed to the F99-template ( Van Essen, 2002 ) and parcellated into 174 regions of interest (ROI) using the LVE atlas ( Lewis and Van Essen, 2000 ). All further analyses were done in R (version 3.6.1; R Core Team, 2019 ).
Nuisance regression was performed using the six motion parameters, their derivatives, their squares, as well as mean white matter and Fig. 1. Competing hypotheses, methods and analytic approaches -a) Alternative accounts of the apparent network fragmentation observed during unconsciousness. The fragmentation account posits a splitting of conscious brain networks into smaller subnetworks during unconsciousness. An alternative muting account explains the apparent fragmentation by a global reduction in correlation magnitude, which results in sparser, more fragmented networks after applying a magnitude threshold. b) We collected twelve five-minute resting state scanstwo at each of six concentrations of isoflurane. The cortex was parcellated using the LVE atlas ( Lewis and Van Essen, 2000 ), and covariance matrices were estimated for each scan. c) We studied dose related changes in network statistics both in the real data, and in surrogate datasets in which either the correlation structure or magnitude (the mean absolute value of the pairwise correlations) were held constant. In the former ( Constant Structure ), we simply scaled the correlation matrix associated with the lowest dose (1% isoflurane) to match the observed mean magnitude for each scan. In the second ( Constant Magnitude ), we scaled all correlation matrices to have mean magnitude equal to the lowest dose. d) To quantify the degree to which the data supported either a constant vs. dose dependent covariance structure, we split each subjects data in two datasets, each comprising one of the two scans at each dose. We then compared the correlation matrices in one half both to the lowest dose ( 1% isoflurane ) or to the corresponding dose ( Dose matched ) in the other.
CSF signals, for a total of 26 regressors. We elected to use the unwhitened residuals rather than performing prewhitening, as the overall residual autocorrelation was minimal (mean Durbin-Watson statistic in the range [1.6-2] for all subjects and all scans; Supplemental Fig. S2).

Functional connectivity estimation and comparison
For each scan, we estimated a covariance matrix from the standardized residuals using the shrinkage estimator proposed by Ledoit and Wolf (2004) , which has been shown to perform better than the sample covariance matrix in high-dimensional, small-sample settings.
Networks were constructed from the real and surrogate datasets by thresholding the entries of each correlation matrix using a one-tailed t-test with a threshold of = . 05 . All graph metrics ( Fig. 3 a) save smallworldness were computed from the binarized correlation matrices using the R package igraph ( Csardi and Nepusz, 2006 ), and community detection was performed using the Louvain clustering algorithm ( Blondel et al., 2008 ). Smallworldness was computed using the R package brainGraph ( Watson, 2019 ).
Structural similarity between correlation matrices was quantified using the distance between the subspaces spanned by the leading five eigenvectors ( Fig. 3 b). These five eigenvectors constitute a basis for a five-dimensional subspace of the full space of 174 ROIs. The set of all such subspaces forms a manifold -called the Grassmann manifold -on which several natural distance measures can be defined ( Edelman et al., 1998 ). We define the distance between two such subspaces to be the arc length of the geodesic between them, equivalent to the magnitude of the vector of principle angles between the subspaces. Specifically, let 1 and

Centering
To center each subject's correlation matrices, we took the approach advocated by Zhao et al. (2018) , which leverages the natural geometry of the space of covariance matrices. We have implemented many of the computations required to replicate the analysis in an R package spdm ( s ymmetric p ositive-d efinite m atrix), which is freely available from a Git repository at https://github.com/areshenk-rpackages/spdm .
The procedure is as follows. For each subject , we computed a geometric mean covariance matrix ̄ using the fixed-point algorithm described by Congedo et al. (2017) , as well as a grand mean ̄ over all subjects and all scans. We then projected each covariance matrix onto the tangent space at the corresponding subject mean ̄ to obtain a tangent vector where log denotes the matrix logarithm. We then transported each tangent vector to the grand mean ̄ using the transport proposed by Zhao et al. (2018) , obtaining a centered tangent vector Finally, we projected each centered tangent vector back onto the space of covariance matrices, to obtain the centered covariance matrix where exp denotes the matrix exponential.
To gain some intuition for this procedure, note that centering a sample of real numbers can be viewed as the process of computing a sample mean, and then applying a translation which takes the sample mean to zero. The resulting centered values can then be viewed as vectors describing the deviations of the observed values from the mean. In the same way, the tangent vectors computed in Eq. (3) can be viewed as the deviations of each correlation matrix from the corresponding subject mean. The transport in Eq. (4) translates these deviations to the grand mean, thus aligning them to a common baseline. A naive approximation to this procedure would involve, for each subject, simply subtracting the mean covariance matrix from the covariance matrix in each scan. This approach has several drawbacks; most importantly, the fact that the difference between two covariance matrices is not necessarily a covariance matrix, and so does not have an obvious interpretation. The procedure we have described above can be viewed as an adaptation of this approach, which respects the geometric structure of the space of covariance matrices.
To visualize the effect of centering ( Fig. 4 a), we derived a twodimensional embedding of the set of covariance matrices before and after centering using uniform manifold approximation (UMAP; McInnes et al., 2018 ). This procedure was applied to the distance matrices computed from the pairwise distances between covariance matrices ( Smith, 2005 ), defined as where ( 1 , … , ) are the generalized eigenvalues of 1 and 2 .

Common component analysis
After centering, we sought an interpretable, low-dimensional summary of the observed covariance matrices in order to characterize changes in functional connectivity across depths of anesthesia. For a single covariance matrix, this could be accomplished by principal component analysis (PCA), where the eigenvectors of the covariance matrix are used as a basis for a low dimensional subspace capturing the dominant patterns of variability in the BOLD signal observed during a single scan. As we had observations for multiple scans and multiple subjects, we considered two approaches for simultaneously decomposing the full set of covariance matrices.
The first is common principal component analysis (CPCA; Flury, 1984;Trendafilov, 2010 ). The CPCA model attempts to simultaneously diagonalize multiple covariance matrices, and so (informally) assumes that the covariance matrices have identical factor structure, though they may differ in the degree to which they express those factors. As this is a restrictive assumption, we also considered a second model -the common component analysis (CCA) proposed by Wang et al. (2011) which relaxes the assumption that the set of covariance matrices may be simultaneously diagonalized at the expense of some interpretability. The CPCA model was fit using the R package cpca ( Ziyatdinov et al., 2014 ), while the CCA model was fit using custom R code implementing the iterative algorithm proposed by Wang et al. (2011) . As both models returned highly similar results, we present only the results of CPCA.

Reduction in correlation magnitude explains apparent network fragmentation
Fig. 2 a displays summary statistics for the correlation matrices estimated from each scan. Consistent with previous work ( Xie et al., 2019 ), increasing dose was associated with an overall decrease in correlation magnitude, with no clear change in the ratio of positive to negative correlations (as in Bettinardi et al., 2015 ). We also observed a decrease in the spectral radius (the magnitude of the largest eigenvalue), suggesting an overall loss of low dimensional structure. This reduction in correlation did not appear to be driven by regions in any single network, but was present both in primary sensory and somatomotor regions, as well as in the frontal and parietal cortices, and connections between them ( Fig. 2 b,c). Although bilateral homologs displayed relatively high functional connectivity compared to non-homologous regions at low doses, this connectivity was likewise seen to sharply decrease with increasing dose ( Fig. 2 d). Fig. 2 e,f displays the decrease in correlation magnitude at the level of individual ROIs, suggesting that all regions show a trend towards decreasing mean connectivity with increasing dose. Together, these effects indicate a global reduction in functional connectivity magnitude as isoflurane dose increases. This effect appears to be associated with an overall decrease in the low-frequency content of the BOLD signal (.01-.03 Hz; Fig. 2 g), which is itself consistent with the observed decrease in correlation, as the autocorrelation introduced by low-frequency signal fluctuations can often result in increased pairwise correlations between time series ( Christova et al., 2011 ). Similar dose-dependent effects were also observed in summary statistics of the raw BOLD signal (Supplemental Figure S1); particularly in the mean and standard deviation, which decrease near monotonically with increasing dose (consistent with Baria et al., 2018;Huang et al., 2014 ), and show an asymptote in the range of 1.75-2% isoflurane, similar to many of the within and between network correlations observed in Fig. 2 a-d.
Next, to determine the extent to which network changes across dose can be attributed to actual changes in network structure, versus an overall decrease in correlation magnitude, we computed graph summary statistics from the thresholded and binarized correlation matrices, as well as from two surrogate datasets in which either the average correlation structure or the correlation magnitude was held constant across dose. These surrogate datasets were important, as they provided a critical basis for interpreting effects in the real data; if the appearance of network fragmentation is due primarily to an overall decrease in correlation magnitude, then the same pattern of results should be observed when

Fig. 2. Correlation magnitude decreases both globally and locally as a function of increasing anesthetic dose -a)
Descriptive statistics for correlation matrices. In all figures, colored lines denote means for each each subject. Solid black line denotes the mean across all subjects. Error bars are one standard error. b) Parcelation of the cortex into primary sensory regions (somatomotor, auditory, and visual), and frontal and parietal cortices. c) Mean absolute correlation within each parcel, as well as between the frontal and parietal cortices (frontoparietal). d) Mean absolute correlation between bilateral homologues and non-homologous ROIs. e) Average muting of functional connectivity per ROI. Figures display the mean absolute correlation, as well as the change in mean absolute correlation relative to the lowest dose (1% isoflurane). f) Surface map of the change in mean absolute correlation at 2.75% isoflurane relative to the lowest dose. g) Mean spectrum at each dose. the exact same correlation structure is held constant across doses, and only the magnitude is allowed to vary. Conversely, fragmentation should be abolished when correlation matrices are scaled to have the same average magnitude. We created the Constant Structure surrogate dataset by replacing each correlation matrix (from 1.00% to 2.75% isoflurane) with a copy of the average of the subject's two correlation matrices at the lowest dose (1.00%), scaled to have matching correlation magnitude to the real data (defined as the mean absolute value over all pairwise correlations). By contrast, we created the Constant Magnitude surrogate dataset by scaling each correlation matrix to have the same average correla-tion magnitude as in the corresponding subject's lowest dose condition (see Fig. 1 c). We then thresholded and binarized each correlation matrix using an uncorrected, one-tailed t-test with a significance threshold of = . 05 . For further comparison, we also binarized the observed correlation matrices using a relative threshold chosen to produce an constant edge density equal to the mean edge density at the lowest dose ( Density Threshold ). If increasing dose is characterized primarily by a overall decrease in correlation strengths, rather than a change in network structure, then graph properties of density thresholded networks should be relatively preserved across dose. Fig. 3. The appearance of network fragmentation is an artifact of an increasingly muted network structure -a) Network summary statistics. Values denote means across all subjects, while error bars denote one standard error. Observed denotes the actual, observed correlation matrices. For Constant Magnitude , correlations at each dose level were scaled to match the average magnitude of the correlations at the lowest dose ( 1% isoflurane). For Constant Structure , the correlation matrix for the lowest dose was replicated across all scans, and scaled to match the observed average magnitude. Networks were constructed by thresholding and binarizing correlation matrices by significance using a one-tailed t-test with = . 05 . For comparison, we also thresholded the observed correlation matrices using a relative threshold to achieve a fixed edge density of .3 ( Density threshold ). Communities were estimated using the Louvain clustering algorithm. Note that the effect of isoflurane on the observed networks is consistent with the scaling of a constant correlation structure. b) Between scan reliability of each network statistic, defined as the correlation between the first and second scans at each dose. Despite the short scan length, network statistics were generally highly reliable. c) Mean spectrum of the correlation matrices at each dose (top). The x-axis is shown on a log scale to better display the effects of dose on the leading eigenvalues. After splitting each subjects data into two halves -comprising the first and second scans at each dose, respectively -we compared the correlation matrices in first half either to the corresponding dose ( Dose Matched ), or to the lowest dose ( 1% Isoflurane ) in the other. Correlation matrices were compared using the distance between the subspaces spanned by the leading five eigenvectors (middle). The bottom figure displays the correlation between the vectorized correlation matrices from the two scans at each dose, suggesting that the reliability of functional connectivity estimates decreases with increases depth of anesthesia.
Broadly consistent with previous work ( Boly et al., 2012b;Hutchison et al., 2014 ) -and the interpretation that brain networks become increasingly fragmented and/or disconnected at increased levels of sedation -we observed an increase in sparsity, the number of communities, modularity, and small-worldness with increasing dose, along with a decrease in network efficiency ( Fig. 3 a). Despite the short scan length, these effects were generally highly reliable, with most subjects showing high correlations in all network statistics across both scans at each dose ( Fig. 3 b).
Notably, we observed a nearly identical pattern of results in our surrogate dataset in which the correlation structure was held constant across dose (Constant Structure), and only the correlation magnitude was varied. Furthermore, these effects were completely abolished in our surrogate dataset in which the correlations were scaled to have common magnitude across dose, or when the networks were thresholded to achieve a fixed edge density. These results, when taken together, suggest that the observed effect of increasing dose on network statistics (i.e., fragmentation) is better explained as a reduction in overall correlation magnitude, or the muting of constant network structure.

Dose effects are well explained by a constant network structure
If the observed dose effects can be explained by the muting of a constant network structure, then the correlation structure at all dose levels should be well approximated by the structure at the lowest dose (1% isoflurane). We directly tested this by splitting each subject's data into two halves, comprising the first and second scans at each dose, respectively. We then compared each whole-brain correlation matrix in the first set either to the corresponding dose in the second, or to the lowest dose, using the distance between the subspaces spanned by the Fig. 4. Latent structure of brain networks is present across all dose levels, but becomes increasingly muted -Common principal component analysis (CPCA) of subjects' centered covariance matrices. The number of displayed components was selected on the examination of the spectrum of the observed correlation matrices (see Figure 3 b). Prior to CPCA, subject covariance matrices were centered to remove static subject differences in functional connectivity. a) Twodimensional embedding of subject covariance matrices by uniform manifold approximation (UMAP; McInnes et al., 2018 ) using a distance matrix constructed by the pairwise geodesic distances. Note the subject level clustering in the uncentered data. b) . Spatial maps for the top five components. c) Component scores for each subject and each scan.
leading eigenvectors (the geodesic distance on the Grassmann manifold; Edelman et al., 1998 ). We chose the subspace spanned by the leading five eigenvectors, as this was the range of the spectrum most strongly affected by dose ( Fig. 3 c, left).
Critically, we found a near identical pattern of results in the lowest dose and in the dose matched comparisons ( Fig. 3 c, middle), indicating that the dominant patterns of network structure at higher dosages are well approximated by the structure already present at the lowest dose.
To visualize this structure, we used common principal component analysis (CPCA; Flury, 1984;Trendafilov, 2010 ) to derive a set of com-ponents summarizing the correlation structure across dosages. Prior to applying CPCA, subject correlation matrices were centered to remove subject differences in functional connectivity. This decision was motivated by previous findings ( Gratton et al., 2018;Xu et al., 2019 ) that variability in functional connectivity in Humans and non-human primates is dominated by stable, subject level effects, which may mask the relatively small differences induced by task manipulations. Consistent with these findings, we found strong clustering of the correlation matrices at the subject level ( Fig. 4 a, left) in the uncentered data, which we were able to remove through our approach ( Fig. 4 a, right). Fig. 4 b,c shows the spatial maps of the components extracted by applying CPCA to the centered correlation matrices, as well as the contributions of these components at each dose. The first of these (CPC1) constitutes a gradient separating visual and medial ventrotemporal areas from parietal, frontal and superior temporal regions. CPC2 constitutes a gradient separating somatomotor areas from the rest of the cortex. CPC3 constitutes a gradient separating right frontoparietal cortex and bilateral superior temporal cortex from the rest of cortex; CPC5 is a near-mirror image of this same gradient; finally CPC4 constitutes a gradient separating bilateral frontoparietal cortex from the rest of the cortex. The top components are consistent with some of the gradients reported by others in both non-human primates ( Margulies et al., 2016;Yacoub et al., 2020 ), and in humans ( Hong et al., 2020 ). Together, these findings support the notion that a constant network structure is present across all dose levels, and that it is only the expression of this constant network structure that changes across dose. With respect to the latter, indeed we find that these components decreased almost monotonically with increasing dose, asymptoting at approximately 2% isoflurane.
Although we are wary of drawing conclusions from a single subject, we also note that subject M5, who experienced an adverse reaction to 2.75% isoflurane, stands out from the remaining subjects both in showing low correlation magnitude compared to the remaining subjects even at low doses ( Fig. 2 ), as well as showing a muted expression of several components ( Fig. 4 ). This may suggest a hypersensitivity to isoflurane which is detectible even at low doses.

Discussion
We analyzed static resting state functional connectivity (rs-FC) under increasing depths of anesthesia in order to characterize dose-related changes in cortical whole-brain network structure. Increasing dose was associated with an increase in network modularity and the number of communities, as well as a decrease in network efficiency, all of which is consistent with the apparent fragmentation discussed in previous literature . However, comparisons with surrogate datasets in which either the correlation structure or magnitude were held constant revealed that these above effects, rather than reflecting a qualitative change in network structure (i.e., network fragmentation), could be fully explained as an overall reduction in correlation magnitude, which we refer to as muting. Further supporting the idea that network structure was unchanged across dose levels, we showed that the principal components of the functional connectivity matrices at each increasing dose level (from 1.00% to 2.75%) were just as similar to the components at the lowest dose (1.00% isoflurane) as they were to the components in another, dose-matched scan. This suggests that there is no explanatory benefit to assuming a change in correlation structure across dose levels, as this structure can be just as well approximated assuming a fixed structure, already present at the lowest dose. Finally, we used common principal component analysis to derive a set of components that summarize the correlation structure across dose, and show that the expression of these components decreases in a near monotonic fashion as dose increases. Taken together, these findings suggest that deepening levels of unconsciousness are associated with the increasingly muted expression of a constant functional network structure, rather than a break-down or fragmentation of this structure. Below, we discuss the implications of these findings; first, to research attempting to characterize whole-brain network changes associated with consciousness using fMRI; and second, to the broader literature attempting to characterize neural signatures of consciousness.

The interpretation of graph theoretic measures in prior work
Several neurobiological theories of consciousness center on the brain's capacity for information exchange across multiple distributed regions, with unconsciousness being a result of the disruption of this information transmission . This perspective has found support from fMRI studies of large-scale brain networks during waking and various states of unconsciousness, which have often relied on the interpretation of graph summary statistics describing the global structure of these networks. This work has found a reduction in network efficiency in both local and global brain networks during anesthetic-and sedativeinduced unconsciousness ( Hashmi et al., 2017;Monti et al., 2013 ), as well as in patients with consciousness disorders, such as unresponsive wakefulness syndrome and the minimally conscious state ( Chennu et al., 2014 ), and has been interpreted as reflecting a dysfunction in communication across the cortex. Likewise, increases in whole-brain modularity and the number of network communities have been observed during both non-rapid eye movement sleep ( Boly et al., 2012b ) and isofluraneinduced anesthesia ( Hutchison et al., 2014;Standage et al., 2019 ), which has been interpreted as reflecting a literal fragmentation of brain networks into smaller, more isolated processing units.
Graph measures are thought to quantify key aspects of information transmission in the brain, and many of them, such as modularity and efficiency, neatly map onto existing frameworks and hypotheses concerning theories of consciousness (and disruptions thereof). However, these measures are two layers of abstraction removed from the underlying functional network structure: First, they are low dimensional summaries of large-scale networks spanning, in many cases, the entirety of the cortex; and second, the networks themselves, often constructed by thresholding BOLD signal correlations, are merely proxies for information flow between cortical regions. The interpretation of these summary measures thus relies on strong assumptions about the ways in which changes in cortical information processing manifest as changes in graph properties. As we have shown, global decreases in functional connectivity, possibly related to changes in the low frequency (.01-.03 Hz) content of the BOLD signal, can manifest as an increase in modularity and in the number of communities in whole-brain networks, despite no apparent selective decoupling between distributed networks (e.g. frontoparietal regions).

Muting and fragmentation are not necessarily mutually exclusive
Our findings demonstrate that global decreases in functional connectivity (muting) can lead to the false appearance of network fragmentation with graph theoretic analyses. Nevertheless, it is important for us to emphasize that such results do not actually undermine the widely held view that unconsciousness is the result of such network fragmentation (e.g. Mashour et al., 2020 ). Indeed, some models argue that consciousness is the result of a non-linear process ( Dehaene et al., 2003;Mashour et al., 2020;Moutard et al., 2015 ), by which there is some minimal level of functional connectivity required in order to sustain effective information transfer between regions. Under this framework, a drop in functional connectivity below the critical threshold required to sustain communication between two brain regions would result in a functional "disconnection " between those regions. In this way, global muting would actually lead to network fragmentation.
The work of Yang et al. (2014) may provide a relevant perspective. The authors study a computational model of resting state functional connectivity and observe a relationship between variability and connectivity, wherein increases in either self-coupling or global coupling strength are associated with increased global signal variability. Huang et al. (2014) likewise found a correlation between BOLD signal variability and BOLD synchronization between local and distant cortical regions. These findings are broadly consistent with our observation that BOLD standard deviation decreases with increasing dose (Supplemental Figure S1), and that the asymptote of this trend in the range of 1.75-2.75% isoflurane matches similar patterns of within and between network correlation strength ( Fig. 2 a-d). The idea of a global decrease in connectivity is thus consistent with both the changes in within and between network correlation we have observed across dose, as well as the changes in global signal variability.
Thus, to clarify, our results are not incompatible with the view that unconsciousness stems from a fragmentation in network structure. What we show, however, is that amidst the backdrop of the global muting of network structure, increases in graph theoretic measures alone, such as whole-brain modularity and the number of network communities, cannot be used as evidence to support this view of network fragmentation.

Local versus long-range connectivity in unconsciousness
One of the main observations to emerge in the anesthesia literature, and which has been frequently used to bolster neurobiological theories of consciousness, is that disruptions in consciousness are associated with a break-down in long-range frontal-parietal connectivity Mashour et al. (2020) . For example, recent work has shown that the decoding of a tactile stimulus is abolished in the macaque primary motor cortex (area M1) but not the primary somatosensory cortex (area S1) after ketamine exposure ( Schroeder et al., 2016 ). This finding was attributed to the preservation of thalamocortical input to S1, but impaired information transfer from S1 to M1. This general result accords with earlier findings showing that connectivity (as measured by Granger causality) between the ventral thalamus and S1, but not M1, is preserved under ketamine Kim et al. (2012) , and it is also broadly consistent with other work (e.g. Ku et al., 2011;Lee et al., 2009 ) showing impaired fronto-parietal connectivity under anesthesia. Here we observed similar reductions in frontal-parietal connectivity under increasing depths of anesthesia; notably, however, we found that this effect was not unique to long-range connectivity, as we also observed comparable decreases within both the frontal and parietal cortices individually, as well as the primary sensory (visual, auditory) and somatomotor (primary sensory and motor) areas.
In light of this seemingly contradictory finding, it is worth noting that the preservation of thalamocortical connectivity under anesthesia may be dependent on the specific anesthetic agent. For example, impaired thalamocortical connectivity has been observed under isoflurane ( White and Alkire, 2003 ), as well as under propofol ( Malekmohammadi et al., 2019 ) anesthesia. Furthermore, disruptions in cortico-cortical connectivity under propofol do not appear to selectively target long-range connections as compared to shortrange connections ( Monti et al., 2013 ). Other researchers have also found disruptions in local neural dynamics under anesthesia, including Vizuete et al. (2014) and  , who find that neural activity in the primary visual cortex becomes more entropic and variable under anesthesia, despite an overall preservation of the repertoire of activity patterns. While these findings, as well as our own results, seemingly depart from the view that long-range frontal-parietal networks play a unique role in conscious experience per se, they nevertheless comport with the broader view, shared by several theories, that unconsciousness stems from a more global disruption of information processing throughout the brain. It is also not inconsistent with views which assign frontal-parietal regions a privileged position due to their role as hub regions facilitating the broadcasting of information widely across the cortex , as this broadcasting may be impaired not only by selective disruption of frontal-parietal regions, but also by more global impairments in connectivity. Finally, it is important to distinguish neural investigations of the effects of anaesthesia with respect to sensory manipulations ( Schroeder et al., 2016 ) versus examinations of spontaneous brain activity at rest. Indeed, whereas Schroeder et al. (2016) observed preserved local neural representations in the presence of a sensory stimulus, examinations of resting-state activity Vizuete et al. (2014) and  reveal disruptions in local spontaneous activity. In this regard, it is worth noting that stimulus onset has been associated with a similar decline in neural variability across various species in both awake and anesthetized conditions ( Churchland et al., 2010 ), and thus the sensory information carried through thalamocortical input may contribute to the robustness of local sensory representations against the local effects of anesthesia.

Methodological considerations
Our findings should be interpreted in light of a few methodological limitations. Before addressing these, we think it is important to stress the distinction between the implications of our findings to theories of consciousness per se, versus the implication of our findings to methodology which attempts to use functional brain networks, estimated by fMRI, to characterize brain-network changes during states of unconsciousness. In the former case, we stress some important aspects of our design that limits our ability to draw general conclusions about the neural signatures of consciousness. In the latter case, we note that the many of our findings are broadly applicable to many similar research designs, and to network-based approaches to the study of consciousness, in particular.
First, although we track network structure across increasing depths of anesthesia, we lack an awake condition for comparison. Our results thus reflect network changes spanning degrees of unconsciousness, rather than the transition between unconsciousness and waking. Bettinardi et al. (2015) observe a gradual re-emergence of correlated brain activity during recovery from deep sedation, culminating in an eventual restoration of waking connectivity. Similarly, Chennu et al. (2014) observe features of functional brain networks which correlate with the degree of wakefulness in patients with disorders of consciousness. These results suggest that our protocol, while lacking an explicit waking comparison, may nonetheless capture a segment of a continuous trajectory spanning wakefulness and deep sedation.
Second, we track network changes only through progressively increasing doses of isoflurane, rather than using a randomized order, or alternating increasing and decreases doses. This decision was partly pragmatic, as we were concerned about residual effects from the extremely high dose in the 2.75% condition (indeed, one subject experienced adverse reactions at this dose), as well as the time required to stabilize subjects when transitioning to a lower dose. Nevertheless, the use of a counterbalanced dose order would allow the study network changes associated with both induction and recovery from anesthesia, although some authors who have tracked brain activity from initial induction to full recovery of consciousness (e.g. Blain-Moraes et al., 2017 ) have observed a gradual return of anesthesia induced network changes back to baseline levels. This suggests that the use of a counterbalanced dose order may not have revealed much new information.
Third, subjects were sedated at 1.5% isoflurane prior to stabilization at 1% before the beginning of the first scan. It is thus possible that residual effects from the initial, higher dose persisted during the initial scans, giving an inaccurate view of brain network structure during lower doses. Although some residual effects are possible, we note that the physiological recordings, BOLD signal correlations, and network summary measures show clear changes in the range of 1-1.5% isoflurane, suggesting that some degree of recovery had indeed occurred, and that network structure at sub-1.5% doses does not merely correspond to residual effects from the initial sedation.
Finally, initial sedation was accomplished through the use of several agents, including propofol, ketamine, and isoflurane, making it difficult to disentangle the effects of isoflurane itself from the distinct effects of the other anesthetic agents. Further complicating interpretation is the fact that different anesthetic agents may affect distinct neurotransmitter systems, and are known to affect whole-brain functional connectivity differently ( Jonckers et al., 2014;Williams et al., 2010 ). First, this kind of polypharmacy is commonplace across numerous studies (e.g. Bettinardi et al., 2015;Hight et al., 2014;Huang et al., 2014;Kim et al., 2012;Li et al., 2019 ), and thus our methodological interpretations remain relevant to the broader literature on this topic. Second, this speaks to a broader issue, which is that research investigating sedative induced unconsciousness faces the hurdle of disentangling brain network changes due to unconsciousness, versus changes due to the unique action of the particular anesthetic compound on the brain. We stress that although this complicates attempts to relate anesthesia-induced changes in brain network structure to consciousness, it does not undermine our broader point that commonly used network summary measures may be entirely misleading in the absence of a detailed characterization of network structure.

Conclusion
Theories relating consciousness to the information processing capacity of the brain must ultimately generate concrete predictions describing the patterns of brain activity supporting conscious awareness. Functional neuroimaging plays a central role in the testing of these predictions, as these techniques allow for the precise characterization of brain activity across various states of consciousness and unconsciousness. Although fMRI affords the opportunity to study the large-scale, structural features of whole-brain networks, the complexity of this data necessitates careful analysis, and it is essential to determine rigorously the way in which predicted changes in the brains functional network structure may manifest as changes in estimated functional connectivity, or in network summary measures. As we have shown, quantitative changes in overall functional connectivity (muting) are sufficient to produce changes in network structure more commonly interpreted as reflecting the selective disconnection of distributed networks. Distinguishing between these possibilities requires the detailed analysis of whole-brain network structure beyond that which is provided by graph summary measures. It likewise requires consideration of the underlying changes in the brains functional circuitry which may give rise to these effects measured by fMRI functional connectivity. Although we characterize the effects which we have observed as muting, we do not address its causal role in consciousness, or whether it may simply mask structural network changes associated with unconsciousness. Recent work ( Pal et al., 2020 ) suggests that anesthesia induced unconsciousness can be dissociated from the effects of anesthesia on the cortex, and thus the relationship between consciousness and cortical information processing may be more complex than can be decoded from BOLD signal correlations alone.