Multi-modal analysis of functional connectivity and cerebral blood flow reveals shared and unique effects of propofol in large-scale brain networks.

Anesthesia-induced changes in functional connectivity and cerebral blow flow (CBF) in large-scale brain networks have emerged as key markers of reduced consciousness. However, studies of functional connectivity disagree on which large-scale networks are altered or preserved during anesthesia, making it difficult to find a consensus amount studies. Additionally, pharmacological alterations in CBF could amplify or occlude changes in connectivity due to the shared variance between CBF and connectivity. Here, we used data-driven connectivity methods and multi-modal imaging to investigate shared and unique neural correlates of reduced consciousness for connectivity in large-scale brain networks. Rs-fMRI and CBF data were collected from the same subjects during an awake and deep sedation condition induced by propofol. We measured whole-brain connectivity using the intrinsic connectivity distribution (ICD), a method not reliant on pre-defined seed regions, networks of interest, or connectivity thresholds. The shared and unique variance between connectivity and CBF were investigated. Finally, to account for shared variance, we present a novel extension to ICD that incorporates cerebral blood flow (CBF) as a scaling factor in the calculation of global connectivity, labeled CBF-adjusted ICD). We observed altered connectivity in multiple large-scale brain networks including the default mode (DMN), salience, visual, and motor networks and reduced CBF in the DMN, frontoparietal network, and thalamus. Regional connectivity and CBF were significantly correlated during both the awake and propofol condition. Nevertheless changes in connectivity and CBF between the awake and deep sedation condition were only significantly correlated in a subsystem of the DMN, suggesting that, while there is significant shared variance between the modalities, changes due to propofol are relatively unique. Similar, but less significant, results were observed in the CBF-adjusted ICD analysis, providing additional evidence that connectivity differences were not fully explained by CBF. In conclusion, these results provide further evidence of alterations in large-scale brain networks are associated with reduced consciousness and suggest that different modalities capture unique aspects of these large scale changes.

Introduction

Consciousness is a hallmark of normal human life, yet its breadth and complexity make consciousness remarkably difficult to study. Anesthetic agents are designed to suppress levels of consciousness and such agents provide a powerful tool to relate changes in consciousness to changes in brain properties (Mashour, 2006). With a few exceptions such as ketamine, anesthetics suppress baseline metabolism (Franks, 2008), suggesting a brain energetic pathway for consciousness (Shulman et al., 2009). Similarly, electroencephalogram (EEG) data shows a transition from high-frequency electrical activity to slow cortical waves as anesthesia takes effect, reflecting a reduction in neuronal firing rate (Alkire et al., 2008). For both metabolism and electrical activity, the greatest changes are observed large-scale brain networks spanning cortical and subcortical regions (Franks, 2008; Hudetz, 2012).

Changes in functional connectivity also point towards similar disruption of large-scale brain networks in anesthesia-based reduced consciousness (Boveroux et al., 2010; Gili et al., 2013; Liang et al., 2012; Liu et al., 2013; Martuzzi et al., 2010; Monti et al., 2013). Cortical regions commonly reported include regions of the default mode network (DMN; Greicius et al., 2008; Martuzzi et al., 2010; Mhuircheartaigh et al., 2010; Stamatakis et al., 2010), the frontoparietal network (FPN: Boveroux et al., 2010; Schrouff et al., 2011), and motor-sensory networks (Martuzzi et al., 2010; Peltier et al., 2005).

However, studies of functional connectivity disagree on which large-scale networks are altered or persevered during anesthesia, making it difficult to find a consensus among studies (Hudetz, 2012). Many functional connectivity techniques applied to investigations of anesthetics have focused on a limited number of seed regions or a predefined network (such as a component derived from independent component analysis), making it difficult to detect changes in more than one network. Focusing on pre-defined specific regions or canonical resting-state networks potentially yields an incomplete picture of the large-scale network changes due to anesthesia-related changes in consciousness.

Additionally, connectivity studies often employ only a single modality (typically resting-state functional magnetic resonance imaging or rs-fMRI) and do not account for changes in metabolism or blood flow. Variations in these factors can modulate the rs-fMRI signal (Jo et al., 2010; Khalili-Mahani et al., 2014; Liu, 2013), leading to mischaracterization of the shared or unique aspect of functional connectivity correlates of reduced consciousness.

In this work, we used data-driven connectivity methods and multi-modal imaging to investigate neural correlates of reduced consciousness for connectivity in large-scale brain networks. Rs-fMRI and CBF data, as measured by pulsed arterial spin labeling (ASL), were collected from the same subjects during an awake and a deep sedation condition. We measured whole-brain connectivity using the intrinsic connectivity distribution (ICD (Scheinost et al., 2012)), a method not reliant on pre-defined seed regions, networks of interest, or connectivity thresholds. The shared and unique variance of connectivity and CBF were investigated using a high-resolution, 268 node functional atlas (Finn et al., 2015). Finally, to account for the shared variance, we present a novel extension to ICD that incorporates cerebral blood flow (CBF) as a scaling factor in the calculation of global connectivity, labeled CBF-adjusted ICD. We hypothesized that unique connectivity differences due to reduced consciousness will span multiple large-scale brain networks when accounting for the shared variance with CBF.

Materials and methods

Subjects

Thirty-two healthy subjects (19–35 years) underwent imaging that included resting state blood oxygenation level dependent (BOLD) fMRI and pulsed arterial spin labeling (PASL) acquisitions during an awake condition and a propofol anesthesia condition within the same imaging session. Fig. 1 shows a timeline of image acquisition. The BOLD fMRI scans were always acquired before the ASL scans. The sample included 19 males and 13 females with an average age of 23.3 ± 3.1 years. Subjects on psychoactive drugs or centrally acting medications and those with a history of renal disease were excluded. MRI experiments were performed under research protocols approved by the Institutional Review Board of Yale University. All patients fasted for more than 8 h before the study.

Fig. 1

Schematic of image acquisition timeline

Baseline resting-state fMRI (rs-fMRI) and cerebral blood flow (CBF) data were acquired while subjects were awake in the resting state. Next, within the same imaging session, a low dose of propofol was administered to the participant with a maintenance infusion target plasma concentration of 2 μg/ml. After allowing the subject to reach an anesthetized state (approximately 10 min), the same sequence of CBF and rs-fMRI data was acquired.

Propofol administration

Propofol was infused to induce an anesthetized state. For propofol administration, an intravenous (IV) cannula was inserted under sterile conditions, and an IV infusion line was used for maintenance infusion (lactated ringer at 100 ml/h). Propofol infusion was administered through a target-controlled infusion (TCI) pump (Stanpump, Stanford University, Palo Alto, CA), which is based on a kinetic model that uses the age, sex, weight and height of the subject as model inputs. The target plasma concentration for TCI was 2 μg/ml plasma level. This concentration produces unresponsiveness to verbal commands and global loss of memory (Purdon et al., 2009) and corresponds to “Deep Sedation” as defined by American Society of Anesthesiologist (ASA) continuum of sedation. As a confirmation of level of consciousness, all subjects did not respond to verbal call. Painful stimuli were not performed to arouse the participants.

For the deep sedation condition, the MRI scans began at least 10 min after the start of propofol infusion. To assess the variability in propofol delivery as determined by plasma level for each subject, two blood samples (5 ml each) were drawn on 23 of the 32 participants. One sample was drawn at least 10 min after the start of propofol infusion and the other right before the end of the propofol infusion. For subject safety, physiological parameters were monitored throughout the imaging session. Blood pressure was monitored every 5 min using a standard blood pressure cuff. Respiratory rate, heart rate, end-tidal CO2, and blood hemoglobin oxygenation were continuously monitored with standard ASA monitors. Oxygen was delivered to the subjects through a nasal cannula at a rate of 4 l per minute for both the awake and deep sedation conditions. During the experiment, subjects were instructed to lie still in the scanner with their eyes closed.

BOLD fMRI acquisition

Magnetic resonance imaging data were acquired on a 3 T whole-body scanner Siemens TIM Trio with a 12-channel phased array head coil. During the awake and the propofol conditions, two functional BOLD runs of 210 volumes each were acquired using a T2*-sensitive gradient-recalled, single-shot echoplanar imaging pulse sequence (TR=2 s, TE=30 ms, FOV=256×256 mm2, flip angle=90°, matrix size 64×64). Each volume consisted of 33 AC-PC aligned slices, with a slice thickness of 4 mm and no gap.

PASL CBF acquisition

The Quantitative Imaging of Perfusion using a Single Subtraction (QUIPSS) pulsed arterial spin labeling (PASL) sequence was used for measuring resting-state CBF in the anesthesia-free and anesthesia conditions (Luh et al., 1999). A slab-selective (100 mm) hypersecant inversion radiofrequency (RF) pulse was used for labeling the in-flow arterial blood water. The RF pulse was applied to a slab 25 mm inferior to the imaging slab for labeling, and the same RF pulse was applied to a slab 25 mm superior to the imaging slab to control for off-resonance effects. The paired labeling and control images were acquired in an alternating manner. A 20-slice ASL acquisition was implemented, and all slices were acquired parallel to the AC-PC line and positioned to provide full brain coverage. The ASL acquisition parameters were: field of view FOV=256×256 mm2; matrix=64×64; bandwidth=2004 Hz/pixel; slice thickness=5 mm with an interslice gap of 2.5 mm. The repetition time was TR=3000 ms; the echo time was TE=26 ms. The acquisition of each slice took approximately 60 ms, therefore the post-labeling inversion time for each slice i, i=1, 2, …, TI(i)=1400+60 (i−1) ms, which was used in CBF quantification, resulting in a post-labeling data acquisition window from 1.4 to 2.6 s. A bipolar gradient of encoding velocity Venc=5 cm/s was applied to the imaging slices to suppress the signal contamination from the labeled arterial water within large vessels. Proton density weighted images were collected using the same perfusion sequence, except for the following changes: TR was set to 10 s; the delay time TD was set to 0 ms; and the inversion time TI was set to maximum to allow a full longitudinal magnetization recovery. Two PASL runs, each consisting of 250 volumes, were acquired for both the awake and the propofol condition.

Anatomical imaging

In addition to the functional data, three anatomical scans were acquired. For multi-subject registration of the fMRI data, 2D T1-weighted images were acquired (TR=300 ms, TE=2.43 ms, FOV=256×256 mm2, matrix size 256×256, flip angle=60°), with the same position and orientation of slices as the functional volumes. For multi-subject registration of the perfusion data, 2D T1-weighted images were acquired (TR=300 ms, TE=3.69 ms, FOV=256×256 mm2, matrix size 256×256, flip angle=60°), with the same position and orientation of slices as the ASL volumes. Finally, during the propofol administration (while waiting for subjects to attain a steady-state propofol condition), a high-resolution 3D Magnetization Prepared Rapid Gradient Echo (MPRAGE) volume was acquired (176 contiguous sagittal slices, voxel size 1 mm3, FOV=256×256 mm2 TR=2530 ms, TE=3.34 ms, flip angle = 7°).

CBF calculations

Perfusion-weighted and the proton-density images were motion corrected using SPM8. Time series of the perfusion-weighted images were obtained by pair-wise “surround” subtraction between interleaved label and control pairs resulting in a temporal resolution of 2TR (Aguirre et al., 2002; Wang et al., 2003; Wong et al., 1997). Perfusion-induced difference maps (ΔM) were calculated by averaging these time series for a condition (either awake or propofol). The mean image of the motion-corrected proton density images ( ) was estimated by averaging multiple acquisitions. The absolute CBF (f) (ml/100 g/min) for this condition is calculated on a per-voxel basis from the following: , where λ is the tissue-blood partition coefficient for water, and accounts for the exchange of labeled magnetization from intravascular to extravascular space and clearance of the labeled blood water out of the capillary bed, τ is the arterial transit time, and T1 is the apparent longitudinal relaxation time. Parameters used in CBF quantification include: T1 = 1490 ms, λ = 0.9 ml/g, α = 0.95, t − τ = 700 ms, and TI=1400 ms for the first slice.

Connectivity preprocessing

The first ten volumes of each functional run were discarded to allow for the magnetization to reach a steady state. Slice-timing and motion correction was done using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/). Images were iteratively smoothed until the smoothness of any image had a full-width half maximum of approximately 6 mm (Friedman et al., 2006; Scheinost et al., 2014) using AFNI’s 3dBlurToFWHM (http://afni.nimh.nih.gov/afni/). This iterative smoothing reduces motion related confounds. All further analyses were performed using BioImage Suite (Joshi et al., 2011). Several covariates of no interest were regressed from the data including linear and quadratic drifts, mean cerebral-spinal-fluid (CSF) signal, mean white-matter signal, and mean gray matter signal. For additional control of possible motion related confounds, a 24-parameter motion model (including six rigid-body motion parameters, six temporal derivatives, and these terms squared) was regressed from the data. The data were temporally smoothed with a zero mean unit variance Gaussian filter (approximate cutoff frequency=0.12 Hz). A gray matter mask was applied to the data, so only voxels in the gray matter were used in further calculations. Finally, for each participant, all preprocessed resting-state runs of the same experimental conditions (awake or deep sedation) were concatenated.

Finally, as group differences in motion have been shown to confound connectivity studies (Van Dijk et al., 2012), we calculated the average frame-to-frame displacement for each participant’s data. In line with current reports, participants with an average frame-to-frame displacement greater than 0.15 mm for any run in either the awake or propofol condition were removed from the analysis (n=11). There were no significant differences for motion between the awake and deep sedation conditions (awake: 0.077 ± 0.04, deep sedation: 0.077 ± 0.03, p=0.98). Finally, we regressed a 24-parameter motion model and used an iterative smoothing algorithm to minimize any motion confounds (Satterthwaite et al., 2013; Scheinost et al., 2014; Yan et al., 2013).

Intrinsic functional connectivity

After preprocessing, functional connectivity of each voxel as measured by the intrinsic connectivity distribution (ICD) was calculated for each subject as described previously (Scheinost et al., 2012). Similar to most voxel-based functional connectivity measures, ICD involves correlating the time series for any voxel with every other time series in the brain or brain hemisphere, and then a summary statistic based on the network theory measure degree was calculated. ICD avoids the need for choosing an arbitrary connectivity threshold by modeling the entire distribution of correlation thresholds using a Weibull distribution: , where r is a correlation between two time series, α is the variance parameter, and β is the shape parameter. This parameterization is akin to modeling the change in network theory metric degree, as the threshold used to calculate degree is increased, with a stretched exponential: , where τ is the correlation threshold, and α and β are the parameters as above.

Specifically, the time series for any gray matter voxel was correlated with every other voxel in the gray matter. A histogram of these correlations was constructed to estimate the distribution of connections to the current voxel. This distribution was converted to a survival function, and the survival function was fitted with a stretched exponential with unknown variance. As variance controls the spread of the distribution of connections, a larger variance indicates a greater number of high correlation connections. Finally, this process is repeated for all voxels in the gray matter resulting in a whole-brain parametric image summarizing the connectivity of each tissue element. ICD was performed on each condition independently resulting in two parametric images of connectivity (one for the awake condition and one for the deep sedation condition) for each participant. The two parametric images for each participant were used for further paired statistical analysis.

CBF-adjusted intrinsic connectivity distribution

Based the results present below, CBF and ICD showed specific regional changes due to anesthesia-based reduced consciousness and these regional changes in CBF and ICD connectivity were significantly correlated. As such, resting state time series in different regions might have different amplitudes and correlation structures across the brain due to changes in CBF. To account for these confounds, we propose a novel weighting scheme to ICD, labeled CBF-ICD. In this analysis, connections between voxels were de-weighted by their CBF, reducing the impact of high CBF on the calculation of whole brain connectivity. Formally, instead of constructing a histogram of correlations as in ICD, a weighted histogram of correlations is used for CBF-adjusted ICD. The weighting factor for the count of each correlation between two voxels is the inverse of the product of the voxel’s CBF ( , where f1 and f2 are the CBF values for the two voxels, and f̄ is the mean CBF value in the gray matter). This multiplicative weighting produces weights that are, experientially, between 0 and 1. Correlations between voxels with lower CBF are weighted near 1 (yielding the same contribution to connectivity as standard ICD), while correlations between voxels with higher CBF are weighted near 0 (yielding little contribution to connectivity). Correlations between voxels with CBF near the gray matter average are weighted by, approximately, or . While it is theoretically possible to produce a weight greater than 1 (i.e. f1*f2CBF in the gray matter was not large enough to observe such low CBF values. After creating the weighted histogram, the same procedures for ICD were followed as described above.

Common space registration

Single subject maps were warped to the MNI standard template for group analysis using a series of linear and non-linear transformations. The functional data (BOLD and CBF) were linearly registered to the 2D anatomical images. The 2D anatomical images were linearly registered to the 3D MPRAGE images. The 3D MPRAGE images were nonlinearly registered to an evolving group average template in an iterative fashion. This algorithm iterates between estimating a local transformation to align individual brains to a group average template and creating a new group average template based on the previous transformations (Scheinost et al., 2016a). All transformation pairs were calculated independently and combined into a single transform warping the participant results into common space. This single transformation allows the single participant images to be transformed to common space with only one transformation, reducing interpolation error.

Group analysis

Connectivity and CBF data were processed for each condition (deep sedation or awake) independently. Paired t-tests were used to compare the resulting parametric images between the awake and deep sedation condition. Imaging results are shown at a cluster-level threshold of p < 0.05 family-wise error (FWE) correction as determined by AFNI’s 3dClustSim program (version 16.0.09) using a cluster-forming threshold of p=0.001, 10,000 iterations, a gray matter mask, and a smoothness estimated using 3dFWHMx with the –ACF option.

To investigate the shared variance between CBF and ICD, we used a high-resolution, 268 node functional atlas (Finn et al., 2015) and averaged CBF and ICD values across subjects in these regions. Shared and unique variances were further investigated using eight canonical large-scale brain networks from (Fig. S1; Finn et al., 2015) to explore spatial variations in the coupling between connectivity and CBF. The functional atlas and corresponding networks were defined in Finn et al. (Finn et al., 2015) using a multigraph K-way spectral clustering algorithm (Shen et al., 2013). Pearson’s correlation was used to assess the similarity between conditions and modalities. Correlations were converted to z-values via the Fisher transform and compared with z-tests.

As an exploratory analysis, we investigated the correlation between changes in ICD or CBF and individual differences in plasma concentration of propofol collected at the end of propofol infusion. For each participant, the average difference between the deep sedation and awake condition for each node in the functional atlas was calculated. These averages were then correlated with plasma concentration of propofol. We used this atlas based approach rather than a voxel-wise approach to improve power given the limited number of participants (n=15) with both low motion fMRI data and plasma concentration of propofol. Results for the exploratory analysis are shown at p < 0.05 uncorrected for multiple comparisons.

Robustness of results to global signal regression

Global signal regression (GSR; performed in this work) remains a controversial preprocessing step as GSR changes the correlation structure between brain regions (Saad et al., 2012) may distort group comparisons (Gotts et al., 2013) and may contain meaningful information (He and Liu, 2012; Scheinost et al., 2016b; Yang et al., 2014). To highlight the robustness of our results with respect to preprocessing choice, we repeated all analyses without removing the global signal.

Results

Physiological measurements

End-tidal CO2 and heart rate were not significantly different from the propofol condition compared to the awake condition. Mean blood pressure (MAP) was reduced significantly by propofol, as shown in Table 1. The plasma levels of propofol on the 23 subjects assayed from the blood samples ranged from 641 to 2170 ng/ml. These plasma levels were positively correlated with the end-tidal CO2 pressure for both the awake and the deep sedation conditions (r=0.73, p=0.006; r=0.88, p=0.002, respectively), and were negatively correlated with the MAP for the deep sedation condition (r=0.48, p=0.02). The plasma levels of propofol were not significantly correlated with the sex, age, weight, or height of the subjects (p=0.47, p=0.61, p=0.88, and p=0.11, respectively).

Table 1

Physiological parameters during awake and propofol conditions.

Vitals (mean ± std)	Awake	Propofol	p-value
Et.CO2, mm Hg	34.2 ± 5.2	33.7 ± 6.7	p=0.84
H.Rate, beats/min	58.9 ± 8.8	61.9 ± 10.5	p=0.83
MAP, mm Hg	86.0 ± 8.5	76.7 ± 8.2	p < 0.001

End-tidal CO2 (Et. CO2). Heart rate (H. Rate). Mean blood pressure (MAP)

Comparison of ICD under anesthesia

Regions of significantly (p < 0.05, corrected) greater ICD compared to the whole brain average are shown for in Fig. 2A for the awake condition, and in Fig. 2B for the deep sedation condition. During the awake condition, greatest ICD was observed in the posterior cingulate cortex, medial prefrontal cortex, anterior cingulate cortex, and insula. In contrast, during the deep sedation condition, greatest ICD was observed in the occipital lobe. Direct comparison of the deep sedation and awake conditions revealed significantly (p < 0.05 corrected, Fig. 2C) increased connectivity in motor and visual networks, including the somatomotor cortex, cuneus, fusiform, and bilateral superior temporal sulcus. Decreased connectivity was observed in the DMN and salience network, including the medial and lateral prefrontal cortex (PFC), bilateral insula, posterior cingulate cortex (PCC), and bilateral supramarginal gyrus. In addition, controlling for participant sex, age, or weight did not change in the results.

Fig. 2

Comparison of ICD under deep sedation

For A) the awake condition and B) the deep sedation condition, regions of significantly (p < 0.05, corrected) greater ICD compared to the whole brain average. C) Direct comparison of the deep sedation and awake conditions revealed significantly (p < 0.05 corrected) increased connectivity in motor and visual networks and decreased connectivity was observed in the DMN and salience network.

Comparison of CBF under anesthesia

Regions of significantly (p < 0.05, corrected) greater CBF compared the to the whole brain average are shown for in Fig. 3A for the awake condition, and in Fig. 3B for the deep sedation condition. During the awake condition, greatest CBF was observed in the posterior cingulate cortex, medial prefrontal cortex, and visual cortex. In contrast, during the deep sedation condition, CBF was qualitatively found to be reduced. Direct comparison of the propofol and awake conditions revealed only significantly (p < 0.05 corrected, Fig. 3C) decreased CBF in the FPN, DMN, and visual networks, including the lateral PFC, PCC, intraparietal sulcus, thalamus, and inferior occipital lobe. In addition, controlling for participant sex, age, or weight did not change in the results.

Fig. 3

Comparison of CBF under deep sedation

For A) the awake condition and B) the deep sedation condition, regions of significantly (p < 0.05, corrected) greater CBF compared to the whole brain average were observed in the cingulate and visual cortex. C) Direct comparison of the deep sedation and awake conditions revealed significantly (p < 0.05 corrected) only decreased CBF in the FPN, DMN, visual networks, and thalamus.

Association between regional CBF and regional connectivity under anesthesia

ICD for the awake condition and the deep sedation condition were significantly correlated (r=0.94, p < 0.001, Fig. 4A), as were the CBF for the awake condition and the deep sedation condition (r=0.90 p < 0.001, Fig. 4B). Regional ICD and CBF were correlated across functional subunits for both the awake condition (r=0.73, p < 0.001, Fig. 4C) and the deep sedation condition (r=0.63, p < 0.001, Fig. 4D), indicating a significant shared variance between connectivity and CBF. Under deep sedation, the Visual I network appeared to deviate the most from the predicted linear fit (Fig. 4D), and removing this network reduced the observed correlation between ICD and CBF (r=0.54, p < 0.001). Furthermore, the shared variance between ICD and CBF was significantly reduced under deep sedation (z=2.16, p=0.03). Finally, changes in ICD and CBF under anesthesia exhibited a modest but significant correlation (r=0.25, P < 0.001, Fig. 4E). Together, these results suggest that a significant shared variance existed between the modalities, but that the changes due to propofol were more unique.

Fig. 4

Correlation between regional ICD and regional CBF

A) ICD values for the awake condition and the deep sedation condition were significantly correlated (r=0.94, P < 0.001) B) CBF values for the awake condition and the deep sedation condition were significantly correlated (r=0.90, P < 0.001). Regional ICD and CBF values were highly correlated across functional subunits for both the C) awake condition (r=0.73, P < 0.001) and D) the deep sedation condition (r=0.63, P < 0.001). However, this coupling is significantly reduced under anesthesia (z=2.16, p=0.03). E) The correlation between change in ICD and CBF under propofol was modest but significant (r=0.25, P < 0.001). Different color circles represent different nodes in the functional atlas clustered by eight large-scale brain networks.

When these correlation analyses were performed independently for each of the canonical network, similar results were observed. For all networks, ICD during the awake and deep sedation condition, CBF during the awake and deep sedation condition, ICD and CBF during the awake condition, and ICD and CBF during the deep sedation condition were significantly (p < 0.001) correlated for each network. However, as shown in Fig. 5, the changes in ICD and in CBF were only significantly (p < 0.01) correlated for frontal-temporal subsystem of the default mode network (Andrews-Hanna, 2012), including the dorsal medial PFC, ventral medial PFC, temporal pole, and temporoparietal juction (as shown in Fig. S1). These results suggest that this subsystem was driving the whole brain correlation and further supporting that the changes due to propofol in the two modalities were unique.

Fig. 5

Correlation between change in ICD and CBF under propofol for eight large-scale brain networks

Correlation between change in ICD and CBF were performed independently for each large-scale brain network. Despite the whole brain showing a modest but significant correlation (Fig. 4E, r=0.25, P < 0.001), only a subsystem of the DMN (1st plot) showed a significant correlation between the change in ICD and CBF. The visual association network (last plot) showed a trend-level correlation. Color circles represent nodes belonging to one of the eight large-scale brain networks. Grey circles represent nodes belonging to the remaining seven brain networks.

Comparison of CBF-adjusted ICD under anesthesia

Regions of significantly (p < 0.05 corrected) difference CBF-adjusted ICD between the awake and propofol condition were similar to the regions identified by the ICD analysis (Fig. 2 and Fig. 6). These regions included increased connectivity in the cuneus and decreased connectivity in bilateral insula, and medial and lateral PFC (Fig. 6). As expected from Fig. 5, regions of the DMN that remained significant were in the core DMN rather than the DMN sub-system. Though fewer significant regions were detected with the CBF-adjusted ICD analysis, all the regions detected in the ICD analysis were detected in the CBF-adjusted ICD analysis when a less stringent significance level was applied (not shown). Together, these results suggest that CBF-adjusted ICD is accounting for additional shared variance across between the two modalities, but changes in CBF do not fully explain changes in connectivity. In addition, co-varying for participant sex, age, and weight did not change in the results.

Fig. 6

Comparison of CBF-adjusted ICD under deep sedation

As ICD and CBF showed significant shared variance (Fig. 4), we repeated our ICD analysis using a novel extension of ICD (called CBF-adjusted ICD) that accounts for CBF. While fewer regions exhibited significant differences in CBF-adjusted ICD, these regions were similar to regions detected in the ICD analysis, suggesting that the connectivity differences are not explained by differences in CBF.

Exploratory analysis of individual differences

As shown in the Fig. 7, several nodes showed strong positive correlations (p < 0.05 uncorrected) between change in connectivity and plasma concentration of propofol. These regions were located in similar regions as the increased ICD under deep sedation (see Fig. 2) such as the temporal sulcus, fusiform gyrus, and motor cortex. No correlations between the change in CBF and plasma concentration of propofol were observed.

Fig. 7

Exploratory correlation between change in connectivity and plasma concentration of propofol at the end of the experiment

Several nodes showed strong positive correlations (p < 0.05 uncorrected) between change in connectivity and plasma concentration of propofol. These regions were located in similar regions as the increased ICD under deep sedation (see Fig. 2) such as the temporal sulcus, fusiform gyrus, and motor cortex. No correlations between the change in CBF and plasma concentration of propofol were observed.

As shown in Figs. S2–S6, the results with and without GSR are largely the same.

Discussion

In a group of participants scanned during an awake and deep sedation condition, we used multimodal imaging to investigate the shared and unique neural correlates of reduced consciousness for functional connectivity (measured by ICD) and CBF in large-scale brain networks. We showed that reduced consciousness was associated alteration in several large-scale networks including reduced ICD in the DMN and salience network, increased ICD in the motor and visual network, and reduced CBF in the DMN, FPN, and visual network. Regional ICD and CBF were significantly correlated for both the awake and deep sedation condition, and the correlation was significantly reduced in the deep sedation condition compared to the awake condition. Changes in ICD were significantly correlated with changes in CBF, suggesting shared variance between the modalities. Finally, we proposed a novel method (CBF-adjusted ICD) to incorporate CBF information in the calculation of connectivity. When accounting for CBF, fewer regions exhibited significant differences in connectivity. However, when using less stringent significance levels, all the regions detected in the ICD analysis were detected in the CBF-adjusted ICD analysis, suggesting that connectivity differences are not fully explained by CBF. Together these results provide further evidence that alterations in large-scale brain networks are associated with reduced consciousness and suggest that different modalities capture unique aspects of these large scale changes.

Reduced consciousness alters ICD and CBF in several large-scale brain networks

During the deep sedation condition, both ICD and CBF were altered in multiple large-scale brain networks with prominent reductions in the DMN and the FPN. These regions are primary targets for anesthesia to reduced consciousness (Hudetz, 2012) with many other studies showing similar reductions of metabolism (Akeju et al., 2014; Fiset et al., 1999; Kaisti et al., 2003; Qiu et al., 2008; Veselis et al., 2004) and connectivity (Deshpande et al., 2010; Greicius et al., 2008; Martuzzi et al., 2010; Mhuircheartaigh et al., 2010) with a variety of anesthetics. Comparable CBF decreases in PCC, medial and lateral frontal regions are observed in the post-ictal period of seizure with a loss of consciousness (Blumenfeld et al., 2004). Likewise, comparable decreases in DMN connectivity are observed in coma (Noirhomme et al., 2010) and sleep (Tagliazucchi et al., 2013; Tagliazucchi and Laufs, 2014). As non-anesthetic based models for disruption of consciousness point towards the same pathways for altering consciousness, a reduction in ICD and CBF in the DMN and FPN may be related to a reduction in consciousness rather than directly related to the administration of anesthesia (Heine et al., 2012).

In addition to the DMN, reduced ICD was observed in bilateral insula, a key node of the salience network (Seeley et al., 2007). The salience network, and particularly the insula, is thought to interact with DMN and FPN to direct attention to either internal or external stimuli (Uddin, 2014). In this context, the insula contributes to consciousness as a node responsible for conscious thought processing, and awareness of affect and somatosensation (Blumenfeld, 2012; Craig, 2009; Uddin, 2014). Additionally, the insula is implicated in interoceptive awareness (Critchley et al., 2004) and integration of interoception and visceromotor responses (Seeley et al., 2012). These reflect essential qualities of consciousness (Blumenfeld, 2012; Sanders et al., 2012) and reports of similar reductions in salience network connectivity during anesthesia are beginning to emerge (Bonhomme et al., 2016; Guldenmund et al., 2013). However, it remains unlikely that disruption of only these systems is sufficient to cause reduced or loss of consciousness (Damasio et al., 2013).

We observed strong increased of ICD in motor and visual regions. Other studies have reported consistent connectivity increases in motor regions (Greicius et al., 2008; Martuzzi et al., 2010) and a perseveration of sensory reactivity under light anesthesia (Davis et al., 2007; Kerssens et al., 2005). In contrast, studies using parcellation and network theory have shown overall decreases in connectivity strength during propofol (Monti et al., 2013; Schröter et al., 2012). In both studies, thresholds were applied to the connectivity which may limit comparisons to our voxel-based, and threshold-free ICD measure. Additionally, parcellation-based analyses only approach voxel-based results when the number of parcels or ROIs is greater than 1,000 (Craddock et al., 2012; Shehzad et al., 2014).

Finally, we observed reduced CBF in the thalamus during the deep sedation condition, consistent with previous studies (Hudetz, 2012). Reductions in consciousness have been hypothesized to occur when the thalamus and cortex become functionally disconnected (Alkire, 2008; Blumenfeld, 2012; Englot et al., 2009).

Shared and unique variance between connectivity and CBF during reduced consciousness

In line with previous reports (Jann et al., 2015; Li et al., 2012; Liang et al., 2013), the spatial patterns of CBF and ICD connectivity are highly correlated suggesting significant shared variance between blood supply and brain functional topology. Under propofol, the shared variance between ICD connectivity and CBF was significantly reduced providing evidence that the relationship between ICD connectivity and CBF is modulated by anesthesia, similar to other pharmacological agents such as caffeine and nicotine (Duncan and Northoff, 2013). Nevertheless, the shared variance between ICD connectivity and CBF is not eliminated under propofol.

Despite this shared variance, changes in ICD were only significantly correlated with changes in CBF for a sub-system of the DMN, suggesting that connectivity and CBF were detecting largely unique, rather than shared, differences due to propofol-induced deep sedation. This result was paralleled with our CBF-adjusted ICD analysis, which produced similar, albeit less significant, results as the ICD analysis. In other words, the gross pattern connectivity differences did not qualitatively change when accounting for CBF.

BOLD fMRI is an indirect measure of neural activity, derived from changes in oxygen metabolism, CBF, and blood volume in response to preceding neural activity (Li et al., 2012). As such, neurovascular factors play a role in resting-state functional connectivity (Jo et al., 2010; Liu, 2013). Pharmacological fMRI studies can alter known neurovascular coupling between spontaneous neural activity and BOLD fluctuations when CBF is significantly changed by study conditions (Khalili-Mahani et al., 2014). Altered the neurovascular coupling could result in measured change in the BOLD fluctuations without a true change in neural activity (Behzadi and Liu, 2006; Cohen et al., 2002). As such, underlying changes in connectivity could be driven by artifactual changes in pharmacological fMRI studies with large changes in CBF. By beginning to account for CBF changes, our connectivity results likely reflect true changes in neural activity associated with anesthesia. This unique variance between connectivity and CBF may be due to nonlinear thresholding effects, where increasing changes in neural activity does not produce an accompanying change in hemodynamics (Sheth et al., 2004).

Methodological considerations

Our CBF-adjusted ICD analysis is designed to account for these potentially confounding CBF effects in connectivity analysis. However, it is not possible to fully decouple the shared variance between CBF and connectivity and assign the variance to a single modality. De-weighting the effects of CBF will ultimately de-weight some true differences in connectivity. As such, it is expected that CBF-adjusted ICD will underestimate of the changes in connectivity even when effects for each modality are relatively unique (as suggesting by the current results).

As most connectivity investigations of anesthesia rely on group differences (Hudetz, 2012), individual differences in consciousness are under explored. We present preliminary data linking individual differences in consciousness to changes in connectivity. Interestingly, regions exhibiting correlation between ICD and plasma concentration of propofol were located mainly in the same regions (temporal sulcus, fusiform gyrus, and motor cortex) that exhibited an increase in ICD under deep sedation. In contrast, regions exhibiting decreased ICD during deep sedate did not seem to be modulated by individual differences in propofol concentration. As the regions correlated with propofol concentration play important role in the processing of sensory information, these differences may be related to differences in cortical reactivity often observed in anesthesia studies of consciousness (MacDonald et al., 2015). However, as these are exploratory results, further replication study is needed.

Strengths and limitations

The strengths of this study include the collection of both CBF and rs-fMRI on the same subjects for an awake and a deep sedation condition, allowing for the investigation of the shared and unique variance between the modalities. We used data-driven connectivity methods. Methods that rely on a priori knowledge are limited in their ability to identify altered connectivity beyond the current hypothesis. As such, ICD and similar measures offer complementary information which can not only be used to confirm previous hypotheses, but also to generate new empirical data for theoretical methods (Cole et al., 2010; Gotts et al., 2012; Scheinost et al., 2015), such as detecting differences that span multiple large-scale networks. Participants’ movement is often a concern during studies of altered consciousness. In this study, we used several methods to minimize confound related to participant head movement, including regressing a 24-parameter motion model and using an iterative smoothing algorithm (Satterthwaite et al., 2013; Scheinost et al., 2014; Yan et al., 2013). Additionally, we applied a strict exclusion criterion for head motion (frame-to-frame displacement greater than 0.15 mm), resulting in no differences in motion between the conditions. As such, it is unlikely motion is significantly contributing to our results. Finally, our results were robust to choice of preprocessing strategy (global signal regression), suggesting that our results are not depended on data processing.

Several limitations also exist. Only a single level of reduced consciousness was assessed. Connectivity and CBF patterns likely differ at difference levels of reduced consciousness. Our results may not generalize to other levels of consciousness or a fully unconscious state. We did not acquire functional data at the onset of the administration of propofol, preventing the study of the dynamics associated with a reduction in consciousness (Lewis et al., 2012). While we accounted for CBF changes in our connectivity analysis, we did not account for individual differences in gray matter density or cortical thickness. Gray matter thickness influence factors such as partial volume and may further alter the functional imaging. However, given the paired design of the experiment, gray matter effects would be similar in both conditions and, likely, cancel out. Lastly, as we acquired two functional imaging modalities, the scanning session during deep sedation was longer than most studies and we did not counterbalance the order of the different modalities. As such, any changes in consciousness over the longer scanning could introduce a systematic difference between the two modalities. However, given the continuous monitoring of the participants by ASA standards and the use of a TCI pump, it would be unlikely for a large shift in consciousness to occur during imaging.

Conclusion

In conclusion, this study used data-driven and multi-modal methods to elucidate unique functional connectivity correlates of reduced consciousness while accounting for shared variance with CBF. During deep sedation, connectivity in multiple large-scale brain networks was altered including decreased connectivity in the DMN and salience network and increased connectivity in the motor, and visual networks, consistent with many prior studies (Hudetz, 2012). Together, changes in connectivity in these networks may represent stable correlates of reduced consciousness that are relatively independent of changes in CBF, and may support a role for the DMN in consciousness. However, given the correlational methods used, further study is needed to pinpoint its exact role.