Static and dynamic functional connectivity analysis of cerebrovascular reactivity: An fMRI study

Abstract Background Cerebrovascular reactivity (CVR) is an important aspect of brain function, and as such it is important to understand relationship between CVR and functional connectivity. Methods This research studied the role of CVR, or the brain's ability to react to vasoactive stimuli on brain functional connectivity by scanning subjects with blood‐oxygenation‐level‐dependent (BOLD) functional magnetic resonance imaging (fMRI) while they periodically inhale room air and a CO 2‐enriched gas mixture. We developed a new metric to measure the effect of CVR on each intrinsic connectivity network (ICN), which contrasts to voxel‐wise CVR. We also studied the changes in whole‐brain connectivity patterns using both static functional network connectivity (sFNC) and dynamic FNC (dFNC). Results We found that network connectivity is generally weaker during vascular dilation, which is supported by previous research. The dFNC analysis revealed that participants did not return to the pre‐CO 2 inhalation state, suggesting that one‐minute periods of room‐air inhalation is not enough for the CO 2 effect to fully dissipate. Conclusions Cerebrovascular reactivity is one tool that the cerebrovascular system uses to ensure the constant, finely‐tuned flow of oxygen to function properly. Understanding the relationship between CVR and brain dynamism can provide unique information about cerebrovascular diseases and general brain function. We observed that CVR has a wide, but consistent relationship to connectivity patterns between functional networks.


| INTRODUC TI ON
Cerebrovascular reactivity (CVR) reflects the brain's ability to modulate cerebral blood flow (CBF) in response to vasoactive stimuli.
CVR reflects dilation and constriction capacity of blood vessels in response to vasoactive stimuli. The dilation and constriction serve to regulate the cerebrovascular blood flow (CBF). Studying CVR provides an important perspective on how the brain functions in conjunction with the vascular system, which could lead to a greater understanding of cerebrovascular disease (Fierstra et al., 2013). CVR can be measured by inducing vasodilation, such as inhalation of a CO 2 gas mixture, while monitoring perfusion-sensitive MRI signals such as blood-oxygenation-level-dependent (BOLD) MRI . The use of gas inhalation to study perturbations in cerebral blood flow as reactions to stimuli dates back to the 1940s (Shenkin, Scheuerman, Spitz, & Groff, 1949) and has been a popular tool since.
The gas mixture used has also changed over time, including both NO 2 -heavy (Kety & Schmidt, 1945, 1948 and CO 2 -heavy (Shenkin, Novak, Goluboff, Soffe, & Bortin, 1953) combinations. Other researchers have studied the relationship between CVR and FNC prior to our research (Tak, Polimeni, Wang, Yan, & Chen, 2015). Cerebrovascular reactivity is typically measured using BOLD signal at the voxel level (Liu, De Vis, & Lu, 2019;Lu et al., 2014) by conducting a linear regression between voxel-wise time course (TC) of the BOLD signal and end-tidal (Et) CO 2 , which is the CO 2 content in the exhaled air and an estimate of arterial CO 2 level in an individual's central nervous system. This method of studying CVR through the use of CO 2 inhalation has been used in previous research to study cerebrovascular disease Yezhuvath et al., 2012) and to study brain networks related to CVR (Liu et al., 2016). It should be noted, however, that there are other techniques to achieve the same goal, such as breath-holding techniques (Bright & Murphy, 2013;Chan, Evans, Rosen, Song, & Kwong, 2015;de Boorder, Hendrikse, & van der Grond, 2004) and venous refocusing for volume estimation (VERVE) (Chen & Pike, 2009;Hoge et al., 1999;Stefanovic & Pike, 2005). VERVE is a newer technique than CO 2 gas inhalation, uses metabolically induced deoxyhemoglobins to deoxygenate subjects' blood and modulate the BOLD signal. However, CO 2 gas inhalation is still a widely used and exhaustively studied technique for measuring CVR.
In this research, we expanded the concept of measuring CVR from the per-voxel level to a per-network level. The per-network CVR provides an understanding of pathophysiology as it relates to functional networks. Capturing the direct relationship between CVR and functional networks could provide deeper insight into how between-network connectivity is altered, moving beyond spatial patterns to provide information about the ongoing dynamics. We accomplished this by calculating the correlation between the EtCO 2 TCs and each network TC and then averaging this correlation across all subjects.
Our research further sought to develop a better understanding of how CVR and vasodilation relate to intrinsic connectivity networks (ICNs) as well as the virtually unstudied area of how vasodilation due to CO 2 inhalation impacts the connectivity relationships between ICNs. We accomplished this by observing and analyzing the influence of the CO 2 inhalation effect on functional network connectivity (FNC), estimated as the cross-correlation between ICN time courses via independent component analysis (ICA). Independent component analysis is an effective tool as it is data-driven and, as such, preserves the vascular relationship within and between networks (Anderson et al., 2011), which is why it was chosen to analyze the CO 2 task data. From the CO -2 task data, we divided the data into CO 2 and room-air time courses based on the EtCO 2 . Static FNC was applied to the CO 2 and room-air segmented ICN time courses in order to capture the overall changes of the FNC time courses, while dynamic FNC was used to capture more nuanced changes that may have been lost while using the full time course to measure FNC. In order to include the time-varying information during the subjects' transitions between the different air-composition intervals, we segmented the dFNC time courses, as opposed to the ICN time courses prior to dFNC. The results of both methods were captured and compared between the room-air and CO 2 conditions. Our dFNC analysis, which was a first for the field, provided nuanced information about the dynamics of the data, which has not been captured by previous research. This analysis of network connectivity, from both static and dynamic perspectives could lead to a better understanding of how the human brain operates in the context of CVR.

| Data collection and acquisition
The study was approved by the Institutional Review Board (IRB) of the University of Texas Southwestern Medical Center at Dallas. Each participant gave written informed consent. The dataset (Hou et al., 2019) included 54 healthy participants that were part of the Dallas Lifespan Brain Study (DLBS)  and consisted of 22 males and 32 female young adults with ages ranging from 20 to 39. During the scans, each subject inhaled a CO 2 gas mixture for approximately 60 s and then spent another 60 s breathing room air.
This cycle recurred for a total of three times (Figure 1).

| Gas delivery system
Each subject was fitted with a nose and mouth apparatus affixed with a non-rebreathing valve, which ensures one-way gas flow, so they could breathe either room air or the CO 2 gas mixture. This two-way breathing valve is affixed to the MRI head coil and connected to a bag containing the CO 2 gas mixture. A research assistant was present in the scanner room to manually operate a control valve to switch between room air and CO 2 gas. A signaling bar was inserted through a wave guide between the control room and the magnet room as a way to communicate with the research assistant when to switch the control valve between CO 2 gas and room air.
The EtCO 2 concentration in the subject's lungs was recorded using a capnograph device. The heart rate and arterial oxygen saturation of the subject were monitored using a pulse oximeter .
During the CVR scan, the subjects began by breathing normal room air for 60 s (30 time points) and then breathed a gas mixture with an elevated level of CO 2 for another minute. From this point forward, the subjects cycled inhalation methods every 60 s (30 time points).
The scan paradigm is shown in Figure 1, which shows the breathing task in detail. CO 2 gas mixture consisted of 21% O 2 , 74% N, and 5% CO 2 . This mixture has an elevated level of CO 2 , compared with the normal room air, which is generally negligible, but with similar levels of O 2 and N 2 (Liu et al., 2016).

| Imaging parameters
The data were collected on a 3-Tesla Philips MRI system with a 32-channel head coil. The BOLD fMRI images were acquired using an echo-planar imaging (EPI) sequence with a repetition time (TR) of 2 s, a field of view (FOV) of 220 × 220 × 150 mm 3 , voxel size of (3.4 mm, 3.4 mm, 3.5 mm), 43 total slices, a 64 × 64 matrix, and a total of 207 volumes. The echo time (TE) was 25 ms, and the flip angle was 80°.

| Data preprocessing
Preprocessing was performed primarily within the SPM software (http://www.fil.ion.ucl.ac.uk/spm/) and custom MATLAB (https ://www.mathw orks.com) scripts. We used the INRIAlign toolbox in SPM to correct for subject head motion. Next, we performed a slice-timing correction using SPM. The data were then warped to the Montreal Neurological Institute (MNI) template and resampled to 3 mm 3 isotropic voxels. Next, spatial smoothing with a 6 mm full width at half-maximum (FWHM = 6 × 6 × 6 mm 3 ). Finally, the TC of each voxel was z-scored (variance normalization).
Subject-specific principal component analysis (PCA) was first used to reduce the subject level TC to 200 principal components. Next, the principal components of individual subjects were temporally concatenated, and the group-level PCA was used to reduce the aggregated components from 200 to 100 along the direction of maximal variance across subjects . The infomax algorithm was applied to maximize spatial independence of the group PCA reduced data resulting in a total of 100 independent components. We selected 100 components based on previous research (Damaraju et al., 2014;Griffanti et al., 2014) which shows that 100 ICs is a sufficient amount to properly capture recognized networks. However, there are certainly other options for how many ICs to select, but we prefer to leave a full exploration of this topic for future work.
The ICA algorithm was repeated 20 times and the most central run was selected to ensure stability (Ma et al., 2011). A group information guided ICA (GIG-ICA) approach from GIFT was used to back-reconstruct the subject-specific spatial maps and TCs from group-level independent components (Du & Fan, 2013). The GIG-ICA approach has been shown to be a more effective artifact removal approach than using single-subject ICA prior to the group ICA analysis (Du et al., 2016) and more sensitive to group differences than a spatiotemporal regression approach (Salman, 2017).

| Post-ICA processing
Intrinsic connectivity networks were identified by assessing the spatial maps and TCs of the independent components. One sample t tests for each spatial map were calculated and then those maps were thresholded to obtain the regions of peak activation clusters for each component. ICNs were selected if the peak activations F I G U R E 1 Mean and standard error across individuals of the paradigm of gas inhalation (top) and the concomitant CO 2 (middle) modulation. As well as the mean and standard error across individuals of the BOLD signal time courses (bottom). The x-axis represents time in seconds Time in seconds covered gray matter and showed minimal overlap with vascular, ventricular, or edge regions . The mean spectral power  was calculated for the TC corresponding to a given component. This information, along with a priori knowledge of ICNs was used to select ICNs. This process was made easier by the high-model-order ICA we chose (i.e., this process was easier due to the 100 available components). These procedures resulted in a total of 42 ICNs ( Figure S1). We then assigned each ICN with a specific functional domain, based on our prior work (Allen et al., 2012). The seven domains consisted of subcortical (SC), auditory (AUD), sensory motor (SM), default mode (DM), attention (ATN), visual (VIS), and cerebellar (CB). Once the ICNs were selected, the corresponding TCs were detrended and then despiked for FNC analysis (Du et al., 2016). As it has been shown to be effective in reducing noise, a low-pass band filter (0.15 Hz) was used to preprocess the ICN TCs prior to computing the FNC .
Post-ICA processing methods, including detrending and despiking were applied as extra steps of cleaning to effectively reduce potential noise and obtain the cleanest possible ICN time courses while not removing important information needed for the ICA framework to properly separate the artifacts and neuronal components within fMRI data. Recent work has studied how processing order changes results in certain capacities (Iraji et al., 2019). However, it is worth noting that there is still room to study the impact of the order of detrending on ICA and artifact removal, which is worth further investigation.

| Data partitioning
The EtCO 2 TCs were used to determine in which time points the subjects were exposed to either the CO 2 gas mixture or room air.
The subject-wise TC was thresholded as either below the average or above the average, defining the room-air and CO 2 conditions used in our experiments. To mitigate noise associated with ambiguous time points, or those steps in which the subjects were transitioning between intervals, we also experimented with groups where scans at the beginning and end of each CO 2 and room-air interval were eliminated. The results from these comparisons showed no statistically significant differences from one another and thus we report only on the results using the above/ below mean approach. However, the time points at the beginning and end of each interval might be ambiguous due to the subjects' transition to or from CO 2 inhalation. In order to evaluate the impact of this, we performed the same experiments while removing the first and last time point of each interval. Our results and conclusions were effectively the same as with the case in which all time points were used. The portioning was done after ICA, but the order with respect to the FNC calculations differed based on whether the FNC was static or dynamic. In the case of sFNC, the subject-wise BOLD TCs were separated before the FNC matrices for both groups were individually calculated. In the case of dFNC, the FNC TCs were portioned into the two groups. The act of portioning the TCs after the dFNC calculations was done so as to not bias the FNC TCs by group.

| Network-wise CVR calculation
As this work focuses primarily on functional network analysis, it was pertinent to quantify the effect of CO 2 at the network level.
This informs us as to which networks are most impacted by vascular reactivity. In order to approximate the network-wise CVR, we calculated the correlation coefficients between each network TC and the EtCO 2 TC for every subject. These coefficients were then averaged across all subjects per network and weighted with the network spatial maps to better visualize network-wise CVR. Because the ICN time courses were z-scored, the standard deviation of the TCs is 1. As such, the correlation values are the same as CVR calculations except for a global scaling value, the standard deviation of the EtCO 2 TCs. Due to this similarity, we can rationally use correlation to represent the effect of CVR on the ICNs. We do this as correlation demonstrates the strength of the similarity between the EtCO 2 and ICN TCs.

| Static functional network connectivity (sFNC)
To compute sFNC, the TCs were segmented into either room-air or CO 2 intervals based on the EtCO 2 average for every subject. Then, the pair-wise correlations between ICN TCs were calculated for each subject, which results in a 42 by 42 symmetric FNC matrix. The columns and rows of the correlation matrix were ordered by the aforementioned domains.

| Dynamic functional network connectivity (dFNC)
In addition, a dynamic FNC analysis was performed on all component TCs, including both room-air and CO 2 time points, which includes sliding-window correlation followed by clustering (Allen et al., 2012). The chosen window size was 30 TR (60 s) in steps of 1 TR, consistent with previous work suggesting this is a good tradeoff between over smoothing and sensitivity to noise (Vergara & Calhoun, 2018). To allow for tapering along the edges, each window was defined as a rectangular window of 30 time points, convolved with Gaussian with a 3 TR full width at half max. We estimated covariance from the regularized inverse covariance matrix (ICOV) using the graphical LASSO framework to reduce noise associated with short time series . In order to impose sparsity, we imposed an L1 norm constraint on the inverse covariance matrix. The log likelihood of unseen data was evaluated to optimize the regularization parameter for each subject in a cross-validation framework. The dFNC TCs, or the correlation matrices per time point for every subject were then segmented based on the EtCO 2 thresholding method.

| Clustering
As has been observed in the past, patterns of network connectivity can reoccur within subjects across time and across subjects. Because of this, we used k-means to cluster the FNC windows in order to minimize the distance between members of a cluster and its cluster centroid (Allen et al., 2012). We used the city-block distance as our measure, due to previous research that suggested city-block was more effective than Euclidean (Aggarwal, Hinneburg, & Keim, 2001). The elbow criterion was used to estimate the model order of five clusters. Initially, we clustered a subset of windows (known as subject exemplars) from every subject corresponding to the windows with maximal variance in correlations between component pairs. The exemplars were obtained by calculating the variance in connectivity across all ICN pairs at each window and selecting windows corresponding to local maxima among this variance TC. From this, we clustered the exemplars and calculated the five centroids. These centroids were then used to initialize a clustering of the entire dataset.

| Statistics
From the sFNC results, we computed average matrices across all subjects for both the CO 2 and room-air results. We also computed a paired t-test per ICN pair, for CO 2 versus room air. The same method of comparison was used for each of the five dFNC states between the CO 2 and room-air matrices. From the dFNC matrices, we also calculated several additional analyses. We computed the transition matrices, or the probability of a subject changing from one state to another, between the five states for both the CO 2 and the room-air results, and then compared the two transition matrices with a paired t-test. The mean dwell time (MDT), or how long a subject stayed in a single state without changing states, and the fraction rate (FR), or how often a given state occurred, were also calculated for both CO 2 and room-air results and then compared via a paired t-test.

| RE SULTS
Building upon previous whole-brain functional connectivity work, we estimated and evaluated 42 ICNs and their corresponding time courses using spatial group ICA. Using these networks, we proposed a solution to calculate the network-wise CVR. We then examined

| Network-wise CVR calculations
We applied the network-wise CVR measurement technique on this data, the results of which can be seen in Figure 2. From these results, we can see the difference between voxel-wise CVR and network-wise CVR. The voxel-wise CVR tends to be more prominent in the gray-matter regions of the brain, as that is where much of the brain's vasculature resides. The network-wise CVR, although it does show similarities to the voxel-wise CVR, there were key differences within certain parts of the brain. There are clear areas of low correlation which can be observed in the network-wise CVR maps. This would appear to be caused by a lack of individual networks in those areas. However, there were also regions that were lower in the voxel-wise maps, possibly due to increased noise in the voxel-wise measurements and the multivariate nature of the network-wise CVR maps. Notably, there were prominent networks with low correlation to the CO 2 effect, including components 79 and 55 in the DMN domain, as well as network 9 in the SC domain. Aside from the differences, the highest and most consistent correlations occurred in networks 18, 19, and 23 in the VIS domain and network 54 in the DMN domains, respectively. There were also several SM networks with relatively high correlation to the CO 2 effect. These network names have been identified in Figure 2.

| sFNC results
The sFNC correlation matrices were calculated separately for CO 2 and room-air time points. The paired t-test results provided a F I G U R E 3 The mean sFNC maps for room-air (left) and CO 2 (right) time points. The black lines separate the FNC maps into the seven domains, labeled as follows: subcortical (SC), auditory (AUD), sensory-motor (SM), default-mode network (DMN), attention (ATN), visual (VIS), and cerebellar (CB). The group differences (using paired t tests) between room-air and CO 2 time points. These values are the FDRcorrected negative log of the p-values multiplied by the sign of the t-statistic. These values have been corrected for multiple comparisons via a false discovery rate (FDR) threshold of 0.05

F I G U R E 4
The mean dynamic FNC maps for both room-air (top row) and CO 2 (bottom row) time points. Each column represents the FNC for each state, from state 1 to state 5. Significant cell-wide differences are visible in States 2, 4, and 5 comparison between the CO 2 and room-air matrices across all subjects. The resulting matrix (Figure 3) was corrected for multiple comparisons using the false discovery rate (FDR), thresholded at 0.05.

| dFNC results
States 2, 4, and 5 showed higher brain connectivity in the room-air time points, as opposed to the CO 2 time points. This is expected since, as the oxygenation of the brain increases (due to venous oxygenation), the BOLD signal becomes less sensitive to oxygenation effects caused by neural activity (Boynton, Engel, Glover, & Heeger, 1996). Both room-air and CO 2 time points showed higher connectivity (compared to other pair-wise correlations) within the VIS domain compared with other domain pairs (Figure 4).
The five dFNC states, computed across all time steps regardless of CO 2 content, were compared with paired t-tests across all subjects for each pair-wise correlation. From this, we saw the greatest differences in states 2, 4, and 5. The results can be seen in Figure 5. From this, we can see large differences between room-air and CO 2 time points within the SM and VIS domains in state 2. The SM domain showed higher correlation overall. We also see smaller differences in the SM domain in state 4. Additionally, we see that state 2 most often occurred in the first portion of the experiment, with very little occurrence during the remainder of the experiment. Figure 6 shows the per-time point occurrence of each state averaged across all subjects. This allows us to visualize changes in the dynamic connectivity which is consistent across individuals. We see that state 2 primarily occurs in the first segment of the experiment, before the subjects have inhaled any amount of the CO 2 gas mixture.
This aligns with our FNC results showing that state 2 had a significant difference between room-air and CO 2 time point.

| Transition matrices
After the dFNC states were calculated, we measured the transition probabilities between states. The transition matrices for both room air and CO 2 show high transition probability within states ( Figure 7) and relatively low transition probabilities between states. FDR-corrected paired t-tests show little difference between states room-air and CO 2 transition probabilities, except for within state 2, which aligns with the other results related to state 2.

F I G U R E 5
Paired t tests for each pair-wise correlation of the dFNC maps for (left to right) states 2, 4, and 5 for room-air and CO 2 . The t tests are the negative log of the p-values, corrected with a false discovery rate (FDR) threshold of 0.05, and multiplied by the sign of the t-statistic. States 1 and 3 were omitted because they had nonsignificant differences F I G U R E 6 The percent occurrence of each state across the entire time series (averaged across all subjects). We see that state 2 occurs most often in the first portion of the experiment, before subjects began inhaling the CO 2

| Mean dwell time results
The mean dwell time (MDT), or the average number of consecutive time points a subject is classified as a given state, was also calculated for both room-air and CO 2 time points. A paired t-test was used to compare the MDT per each state between the room-air and CO 2 time points. MDT across all subjects showed significant differences between room air and CO 2 for states 1 and 2 ( Figure 8).

| Fraction rate results
The fraction rate (FR), or the total number of time points a subject is classified as a given state, were also calculated. These results can be seen in Figure 8. As with the MDT results, paired t-tests were used to compare the FR of each state between room-air and CO 2 time points. Similar to the MDT results, the FR of states 1, 2, and 4 were significantly different between CO 2 and room-air time points.

| D ISCUSS I ON
Cerebrovascular reactivity is a powerful approach to study the human brain. In some recent studies, it has been shown that large-scale resting networks can be estimated from CVR data (Hou et al., 2019;Liu et al., 2016). However, CVR data have not yet been studied in the context of dynamic network connectivity and associated dynamics from ICA. Our analysis quantified the effect of CO 2 on both sFNC and dFNC using ICA. This provides an evaluation of the connectivity across the entire brain during a CVR experiment, providing results similar to previous research (Clarisse, Mazerolle, & Jean Chen, 2015;Tak et al., 2015).
We showed, primarily from the dFNC results, that inhaling CO 2 reduces the overall functional network connectivity. These results reinforce and extend previous research (Madjar et al., 2012;. Reduced functional connectivity can be seen both between domain and within domain.
We present novel results showing an estimation of CVR at the network level. This was accomplished by calculating the correlation between each network and EtCO 2 TC, an estimation of CVR.
The EtCO 2 TC is used as an approximation of the CO 2 content in arterial blood within the brain , which acts as the causal agent of vascular reactivity. Currently, there is much discussion to whether the stimuli during hypercapnia are isometabolic or not. Our primary assumption is hypercapnia effects do induce changes in neural activity. Previous work has challenged this suggesting that hypercapnic stimuli are isometric, or that neural activity does not change with respect to the baseline during a F I G U R E 7 The transition matrices between all 5 states for room-air (left) and CO 2 (middle) time points. The paired t tests between CO 2 and room-air time points (right) are the negative log of the p-values, corrected with an FDR correction with a threshold of 0.05. Only one transition cell, 2-2 passed the FDR threshold F I G U R E 8 (left) The FR of all 5 states comparing CO 2 and room-air time points. (right) Mean standard error of dwell times for all five states comparing CO 2 and room-air time points hypercapnia challenge (Chen & Pike, 2009). However, based on other research , we suggest that although the cause and effect is debatable, this does not invalidate the results from CO 2 manipulation. As such, we suggest that high correlation between the EtCO 2 and network BOLD time courses may imply that vascular reserve in a given network is abundant.
Results showed high correlation between the EtCO 2 TC and BOLD signal in both the VIS and SM domains. A benefit of looking at the network-wise CVR relationship is that we can see the interplay between CVR and brain function. This interplay could be explored more in-depth in the future with a more nuanced analysis of the correlation between functional domains and CVR estimation.
Future research could also focus on new techniques to evaluation the relationship between CVR and functional networks that are more nuanced than correlation alone.
Compared with the sFNC results, the dFNC results showed much larger differences between the room-air and CO 2 results, both in local patterns and across the whole brain. This dynamic analysis provides nuanced information about how vasodilation impacts brain connectivity that is not detected in a static analysis. For example, the dFNC analysis captured information about how the brains changed from the beginning of the experiment (before CO 2 inhalation) and the time points during the experiment. The dFNC analysis also appears to be more sensitive to the CO 2 -induced brain activity changes, capturing multiple significant changes in more than one state.
From the dFNC maps, as well as the state occupancy rates (Figure 7), we consider state 2 to be a state in which the room-air time points are mostly free from the CO 2 effect. We observed that state 2 has the most dominant anticorrelation patterns between functional domains, including patterns between the SM and DM domains, the SM and ATN domains, and between the VIS and ATN domains; patterns that have been seen in previous research (Fox et al., 2005, Fox, Zhang, Snyder, & Raichle, 2009Uddin, Clare Kelly, Biswal, Xavier Castellanos, & Milham, 2009). These anticorrelation patterns are most dominant in state 2. Within this state, the CO 2 effect had the highest impact on lowering the connectivity within the sensory-motor domain (seen in Figure 5). This result is consistent with previous findings showing that CO 2 impacts the SM domain (Golestani, Kwinta, Strother, Khatamian, & Jean Chena, 2017;Liu et al., 2013;Mazerolle, Ma, Sinclair, & Pike, 2016), as well as the VIS domain, meaning that our findings that CO 2 inhalation reduces connectivity within these two domains accompanies an overall reduction of activity within the two domains. These results are congruent with results from the network-wise CVR estimation as there was both a significant impact of CO 2 on the FNC of SM and VIS domains, and both domains were highly correlated with the EtCO 2 TCs. We suggest that this congruency adds robustness to our conclusions about both FNC differences and the effectiveness of our network-wise CVR estimation.
Based on our findings that there are marked differences between state 2 and the other states, we show that generally, the BOLD signal after the start of the CO 2 inhalation was affected by the CO 2 inhalation, even during the room-air intervals. This suggests to us that the minute-long periods of room-air inhalation may not be enough time for the average subject to recover from the neural modulation effects of CO 2 inhalation.
State 4 also produced interesting results, in that it showed more global differences than the other states ( Figure 6). This indicates that although CO 2 has significant impact on specific regions, it also has a significant global effect as well. This result may indicate that the vasodilation caused by CO 2 occurs indiscriminately across the entire brain. State 5 shows a similarly global difference, but with a smaller effect.
The MDT and FR are secondary metrics used to evaluate time-varying information of the FNC patterns, which gives us a broader perspective than what FNC maps alone provide. From our results, the FR and MDT show significant differences within state 2 (room air > CO 2 ), which is congruent with the FNC maps, but they also show significant differences within state 1 (CO 2 < room air).
From Figure 6, we speculate that this might be related to the fact that state 1 contains few time points from the initial portion of the experiment. The dFNC matrices show insignificant cell-wise differences within state 1, which implies that both room-air and CO 2 intervals are impacted by the CO 2 effect, as most of the time points occur after or during CO 2 inhalation. However, the FR and MDT do show significant differences, which may be due to the lack of roomair time points clustered as state 1, meaning there are fewer roomair time points compared with CO 2 time points. This could possibly increase the difference within the MDT between room-air and CO 2 time points and would definitely increase the FR differences. We also see that there are significant differences within the state 5 FNC maps between room air and CO 2 , but no significant difference in the FR or MDT. It is also possible that this is, in part, due to the opposite effect found in state 1. Approximately 1% of all time points within the first interval are clustered as state 5. This may be the cause of some of the differences between the FNC maps but may also contribute to the similarities found in the FR and MDT. The first interval is a slightly larger length of time than the other intervals, 80 s compared with 60 s. But, due to the lower number of state 5 time points within the first interval compared to the other intervals, there would be a more equal number of room-air time points and CO 2 time points. This would show the opposite effect from state 1, as it may reduce the MDT differences, and would most likely reduce the FR differences.

| CON CLUS ION
Our network-wise CVR calculation is a simple method to depict the relationship between CVR and individual ICNs. Whether this relationship is causal one way or another is still up for discussion and future research. We suggest that this method could be used in future analyses of CVR to capture more of the spatial relationships using the multivariate connectivity networks. The network-wise CVR maps showed the relationship between CVR and networks, which, when compared to per-voxel CVR maps, show distinctive differences. We showed high network-wise CVR in both the SM and VIS domains.
The results from our experiments showed widespread differences between room-air and CO 2 time points. We saw that across the whole brain, CO 2 time points showed lower network correlation values than the room-air time points. The dFNC analysis, likely a more natural way to analyze brain connectivity, shows more sensitivity to CO 2 effects not detected by the sFNC. For instance, from the dFNC results, we concluded that state 2 was most prevalent prior to exposure to CO 2 . Given this, it may be useful to utilize this connectivity pattern to predict breathing normally or breathing a CO 2 heavy air mixture. Due to the differences between the states, we also suggest that the minute-long period of time between CO 2 inhalation intervals was not enough time for the subjects to fully recover from the CO 2 gas-induced neural modulation. We also concluded that during the task portion of the experiment, the network correlations across the whole brain for both global and local effects, with local effects primarily affecting the SM and VIS domains. The observed effects of CO 2 on the SM and VIS domains are comparable with the network-wise CVR calculations which showed high correlation between BOLD signals in these domains and the EtCO 2 TCs respectively.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.