Investigating intrinsic connectivity networks using simultaneous BOLD and CBF measurements

When the sensory cortex is stimulated and directly receiving afferent input, modulations can also be observed in the activity of other brain regions comprising spatially distributed, yet intrinsically connected networks, suggesting that these networks support brain function during task performance. Such networks can exhibit subtle or unpredictable task responses which can pass undetected by conventional general linear modelling (GLM). Additionally, the metabolic demand of these networks in response to stimulation remains incompletely understood. Here, we recorded concurrent BOLD and CBF measurements during median nerve stimulation (MNS) and compared GLM analysis with independent component analysis (ICA) for identifying the spatial, temporal and metabolic properties of responses in the primary sensorimotor cortex (S1/M1), and in the default mode (DMN) and fronto-parietal (FPN) networks. Excellent spatial and temporal agreement was observed between the positive BOLD and CBF responses to MNS detected by GLM and ICA in contralateral S1/M1. Values of the change in cerebral metabolic rate of oxygen consumption (Δ%CMRO2) and the Δ%CMRO2/Δ%CBF coupling ratio were highly comparable when using either GLM analysis or ICA to extract the contralateral S1/M1 responses, validating the use of ICA for estimating changes in CMRO2. ICA identified DMN and FPN network activity that was not detected by GLM analysis. Using ICA, spatially coincident increases/decreases in both BOLD and CBF signals to MNS were found in the FPN/DMN respectively. Calculation of CMRO2 changes in these networks during MNS showed that the Δ%CMRO2/Δ%CBF ratio is comparable between the FPN and S1/M1 but is larger in the DMN than in the FPN, assuming an equal value of the parameter M in the DMN, FPN and S1/M1. This work suggests that metabolism-flow coupling may differ between these two fundamental brain networks, which could originate from differences between task-positive and task-negative fMRI responses, but might also be due to intrinsic differences between the two networks.

A large body of literature has demonstrated that spontaneous, low frequency fluctuations in BOLD (blood oxygenation level dependent) functional magnetic resonance imaging (fMRI) signals are highly correlated between brain regions that share a common functional specialization (Biswal et al., 1995;De Luca et al., 2006;Fox et al., 2005;Greicius et al., 2003;Lowe et al., 1998;Smith et al., 2009) (see (Cole, 2010;Van Dijk et al., 2009) for recent methodological reviews). The spatiotemporal coherence of BOLD signals is commonly called "functional connectivity" and is used to define intrinsic connectivity networks (ICNs) during both rest (Damoiseaux et al., 2006) and task-driven paradigms (Sadaghiani et al., 2010;Smith et al., 2009).
The importance of brain network dynamics in supporting cognitive function is becoming increasingly apparent (Laird et al., 2011;Leech et al., 2011;Mennes et al., 2011;Menon and Uddin, 2010;Seeley et al., 2007;Spreng et al., 2010). BOLD fMRI provides a means to study how the cooperation of macroscale brain units gives rise to complex behaviour. Furthermore, ICNs also provide a valuable opportunity to study how the functional architecture of the brain changes in disease pathology. In comparison with healthy controls, alterations in resting-state functional connectivity have been reported during normal ageing (Andrews-Hanna et al., 2007) and in a wide range of pathologies, including chronic pain (Baliki et al., 2008(Baliki et al., , 2011, depression (Veer et al., 2010), schizophrenia (Lynall et al., 2010), epilepsy (Bettus et al., 2009;Waites et al., 2006) and Alzheimer's disease (Agosta et al., 2012). The neuro-scientific information extracted in these studies has the potential to be clinically valuable in aiding diagnosis and monitoring the progression of neurological diseases.
The default mode network (DMN)  is the most widely studied ICN. The DMN is preferentially engaged during rest and displays task-induced reductions in BOLD signal in response to the majority of tasks (Hutchinson et al., 1999). These DMN deactivations are commonly observed concurrently with increased BOLD signal in the task-positive or fronto-parietal network (FPN) (Dosenbach et al., 2007;Fox et al., 2005;Spreng, 2012;Vincent et al., 2008). Although consistently identified during the resting-state (Damoiseaux et al., 2006) and widely implicated in task performance (Corbetta and Shulman, 2002;Raichle and Snyder, 2007;Spreng and Schacter, 2012;Spreng et al., 2010), these ICNs are not robustly detected by conventional general linear model (GLM) analysis in all experimental paradigms. This is possibly because the particular attentional demands and level of cognitive engagement (Pallesen et al., 2009) that characterises the putative function of the DMN and FPN can result in subtle responses that display a high degree of trial-by-trial variability and/or non-canonical haemodynamic response shape. However, independent component analysis (ICA) allows coherent spatiotemporal patterns of brain activity to be resolved without a-priori assumptions concerning the timing or shape of the haemodynamic response (McKeown et al., 1998). Obtaining a better understanding of how the activity of the DMN and FPN supports the functional architecture of the brain requires more fundamental characterization of the neurophysiological origins of ICN activity.
Evidence of a neuronal origin underlying ICN activity as defined by BOLD-sensitive fMRI is provided by patterns of coherent fluctuations in neuronal activity during intra-cortical recordings in primates (Leopold and Logothetis, 2003;Scholvinck et al., 2010;Vincent et al., 2007) and also between regions of the DMN in human epilepsy (Jerbi et al., 2010;Miller et al., 2009). These findings have recently been complemented with non-invasive studies of ICNs in humans using magnetoencephalography (MEG) (Brookes et al., 2011). Additionally, functional connectivity is localised to grey matter (De Luca et al., 2006) and often reflects the brain's structural connectivity (Greicius et al., 2009;van den Heuvel et al., 2009). Although MEG (Brookes et al., 2011) and positron-emission tomography (PET)  have been used to measure ICN activity in humans, to-date the majority of studies have used BOLD fMRI. Whilst an understanding of the neuronal origins and functional significance of ICN activity is being progressively established, the complex neurophysiological origin of the BOLD signal means that changes in ICN BOLD activity should be carefully interpreted (Power et al., 2012;Van Dijk et al., 2012). Previous work has improved the understanding of the effect of physiological variability on BOLD signals (Birn, 2012;Wise et al., 2004). This issue is especially important when comparing BOLD measurements of activity between healthy and ageing or diseased brains, where alterations in neurovascular coupling are expected (D'Esposito et al., 2003;Kannurpatti et al., 2010). It is therefore imperative to obtain a more detailed understanding of the underlying neurophysiology and metabolic demand of ICN activity as measured by BOLD fMRI.
A powerful, non-invasive approach for investigating the neurovascular coupling in ICNs is to concurrently record the BOLD signal with arterial spin-labeling (ASL) MRI, which measures the local cerebral blood flow (CBF) or perfusion (Detre et al., 1994). CBF measurements provide a method to probe brain physiology and, when combined with the BOLD signal, allow the estimation of the cerebral metabolic rate of oxygen consumption (CMRO 2 ), a fundamental parameter of the energy supply supporting brain function. Coherent fluctuations in resting-state CBF measurements have already been reported between regions of the visual, motor, attention and default mode networks with good spatial agreement with BOLD data (Chuang et al., 2008;De Luca et al., 2006;Li et al., 2012;Liang et al., 2013;Viviani et al., 2011;Wu et al., 2009;Zou et al., 2009). Taken together, these studies indicate that spontaneous fluctuations in the BOLD signal at rest likely originate from coupled changes in cerebral blood flow and metabolism which reflect changes in neuronal activity. However, although the relationship between changes in BOLD, CBF and CMRO 2 has been investigated in the primary visual and sensorimotor cortical areas that are directly driven by afferent stimulus input (Ances et al., 2008;Chiarelli et al., 2007a;Kastrup et al., 2002;Leontiev et al., 2007;Stefanovic et al., 2004), relatively few studies have investigated these relationships in ICNs beyond these primary sensory regions (Liang et al., 2013;Lin et al., 2010). This is despite widespread evidence of the involvement of ICNs in supporting brain function and task performance (Laird et al., 2011;Meier et al., 2012;Mennes et al., 2011;Sadaghiani et al., 2010;Sala-Llonch et al., 2012).
In the present study, we used concurrent BOLD and CBF measurements to perform a detailed comparison of CMRO 2 changes induced in sensorimotor regions (S1/M1), DMN and FPN network activity during median nerve stimulation (MNS), providing a more complete characterisation of the complex haemodynamic changes that accompany neural activity in these ICNs. We compare the ability of conventional GLM analysis and ICA to identify these ICNs in both the BOLD and CBF data. Furthermore, we compare the spatial localisation and temporal correlation of the extracted responses to MNS between BOLD and CBF measurements, for each ICN. We validate the use of ICA timecourses to calculate the change in CMRO 2 induced by MNS, by comparison with S1/M1 responses calculated using GLM methods. Subsequently, we use ICA to compare the coupling ratio between the percent signal change (Δ%) in CBF and CMRO 2 in the DMN and FPN, and relate this to the Δ%CMRO 2 /Δ%CBF coupling ratio in S1/M1 that is directly driven by afferent MNS input.

Stimulation paradigm and data acquisition
Data were collected on 18 right-handed subjects (age 27 ± 3 years, 8 female). Written, informed consent was obtained from all participants and this study was conducted with approval from the local ethics committee. These data were acquired as part of an EEG-fMRI study, which used different analysis methods to investigate the fMRI post-stimulus undershoot (Mullinger et al., 2013).
Stimulation was applied to the median nerve of the right wrist using square wave electrical pulses of 0.5 ms duration (Digitimer DS7A, Letchworth Garden City, UK). The stimulation current amplitude was set just above the individual motor threshold (range 2.6-7 mA, mean 4.6 ± 1 mA) so as to cause a small thumb distension. Individual blocks of stimulation consisted of 10 s of 2 Hz MNS, followed by approximately 20 s periods of passive rest. The onset of each period of MNS was jittered by up to 500 ms to reduce the expectation of stimulus onset. A total of 40 blocks were acquired in each subject.
A Philips Achieva 3 T MRI scanner (Philips Medical Systems, Best, Netherlands) with an 8-channel receive coil was used to acquire BOLD and ASL data. Cardiac and respiratory cycles were simultaneously recorded using the scanner's vector cardiogram and respiratory belt. S1/M1 was localised during an initial experiment where BOLD data (TE = 40 ms, 64 × 64 matrix, 3.25 × 3.25 mm 2 resolution, flip angle = 85°, SENSE factor = 2, TR = 2 s, 20 slices of 5 mm thickness) were acquired during ten blocks of 2 Hz MNS. IViewBOLD (Philips GLM analysis software) was used to provide real time statistical maps of the BOLD response to MNS in S1/M1. For the main simultaneous BOLD-ASL experiment, ten contiguous axial slices, providing approximately half brain coverage, were centred on the localiser slice exhibiting peak activation in contralateral S1/M1. A FAIR Double Acquisition Background Suppression (DABS) sequence (Wesolowski et al., 2009) was used for simultaneous acquisition of BOLD and ASL data (TE = 13 ms [ASL], 33 ms [BOLD], 2.65 × 2.65 × 5 mm 3 voxels, 212 mm FOV, SENSE factor = 2; TR = 2.6 s for each label or control acquisition, ASL label delay = 1400 ms, background suppression pulses at T BGS1 /T BGS2 = 340 ms/560 ms, scan duration~21 min). To facilitate co-registration and normalisation of functional data, a single EPI volume was acquired with the same geometry as the FAIR DABS imaging (TE = 40 ms, flip angle = 85°, SENSE factor = 2, TR = 8 s), along with both local and whole-head anatomical images, each at 1 mm isotropic spatial resolution.

Data pre-processing
The DABS dataset was separated into BOLD and ASL data for subsequent analysis. Physiological noise associated with the cardiac and respiratory cycles were regressed out of the BOLD data using RETROICOR (Glover et al., 2000), whereas the use of background suppression reduces physiological noise in the ASL data. All data were then motion corrected using FLIRT (Jenkinson et al., 2002) (FSL, http:// www.fmrib.ox.ac.uk/fsl/). Three subjects were excluded from further analysis due to multiple, gross head movements (N3 mm). Data were then interpolated (interp.m function in MATLAB) to an effective TR of 2.6 s, taking account of the timing difference between the ASL and BOLD timecourses (1.2 s, TR-TI). Label-control image pairs of ASL data were subtracted to create CBF images. The BOLD-weighted images of each label-control pair were averaged to produce mean BOLD-weighted data. Further pre-processing was carried out in SPM8 (http://www.fil. ion.ucl.ac.uk/spm/). BOLD and CBF data were smoothed with a Gaussian kernel (5 mm FWHM) and then normalised to the MNI standard brain using affine matrix transforms calculated from the BOLD data.
Functional data analyses were then performed to identify brain networks where the BOLD and CBF signals were modulated by the MNS, first using a conventional univariate GLM to identify voxel-wise, linear correlations between the fMRI signal and the timings of the MNS; and secondly using multivariate ICA to identify coherent spatial patterns of signal fluctuations without a priori assumptions on the timing or shape of the haemodynamic response (HR).

Spatial localisation and response timecourse extraction: GLM
GLM analyses were performed using FEAT 5.98 (www.fmrib.ox.ac. uk/fsl/). At the first-level, individual subject's BOLD and CBF responses to MNS were modelled using a constant-amplitude boxcar regressor of the stimulation timecourse convolved with the canonical doublegamma HRF. Second level, fixed-effects analyses were then performed to calculate group average activation maps for both positive and negative contrasts. BOLD Z-statistic images were threshold using clusters determined by a Z N 2.3 and cluster-corrected significance threshold of p b 0.05. Due to the intrinsically lower contrast-to-noise ratio (CNR) of ASL data, CBF Z-statistic images were less stringently threshold at Z N 1.6, uncorrected. The group-level conjunction of significant positive BOLD and CBF responses to MNS in S1/M1 was calculated and a mask created. The group conjunction mask was then used to mask the individual subject BOLD Z-statistic images and cubic (3 × 3 × 3 voxels) regions of interest (ROI) were centred on the peak voxel. To facilitate subsequent comparison of BOLD and CBF responses, these ROIs were used to extract the mean timecourses for both BOLD and CBF data, for each subject. BOLD and CBF timecourses were epoched based on MNS timings to extract single block (0-30 s) responses. Responses were then converted to percentage signal change relative to the mean of the final 6 s of the block (defined as baseline), and averaged across blocks. The peak BOLD signal change occurring in the first 20 s of the block (the period of MNS-induced signal change) was then found. The mean BOLD response amplitude (Δ%BOLD) was calculated as the average signal within a 2.6 s (1 TR period) window centred on this peak latency to take into account the interpolation. To allow for slight differences in time-to-peak of the BOLD and CBF signal, whilst ensuring that the equivalent signal change was measured, the peak signal change in the mean CBF response within a time window of ±TR from the timeto-peak of the BOLD response was identified. Finally, the mean CBF response amplitude (Δ%CBF) was calculated as the average signal within a 2.6 s window centred on the peak response latency.

Spatial localisation and response timecourse extraction: ICA
Separately, BOLD and CBF datasets were temporally concatenated across all subjects and MELODIC (Beckmann and Smith, 2004) used to decompose the group data into 20 independent spatial maps and their associated time-courses. Dual-regression (Beckmann, 2009) was then used to identify individual subject timecourses and spatial maps for each independent component (IC). Separately for BOLD and CBF, a single component representing the stimulus response in contralateral S1/M1, as well as the DMN and FPN networks were found. The DMN and FPN were identified from the characteristic spatial pattern of these networks (DMN: precuneus/posterior cingulate (PCC), bilateral inferior parietal lobe (IPL) and medial prefrontal cortex (mPFC); FPN: bilateral intraparietal sulcus (IPS) and lateral prefrontal cortex (lPFC)) (Damoiseaux et al., 2006;De Luca et al., 2006). The following analyses were then performed for contralateral S1/M1, DMN and FPN. Firstly, to assess the spatial similarity of ICNs defined from BOLD and CBF measurements, IC maps were threshold at a Z-statistic N 4 (BOLD) and Z N 2 (CBF), and the conjunction between BOLD and CBF regions calculated at both the group and individual level. Secondly, the correlation between the BOLD and CBF time-series for each subject was calculated as a measure of the temporal similarity between the two signals. As described for the GLM analysis, for each subject both BOLD and CBF timecourses were epoched to extract a single block (0-30 s) response. Responses were then converted to percentage signal and averaged across blocks. Separately for each of the S1/M1, DMN and FPN ICA responses the maximal signal change (whether positive or negative) occurring in the first 20 s of the average block-response, the time-period within which all signal changes peaked on average, was calculated from the mean BOLD response timecourse (Δ%BOLD). The equivalent maximal signal change in the mean CBF response (Δ%CBF) within a time window of ± TR from the time-to-peak of the BOLD response was then found. The mean BOLD and CBF response amplitudes were then calculated as the average signal within a 2.6 s window centred on the latency of these peak responses.
A normalised measure of the single-trial amplitude variability in each network was then obtained to allow comparison of the response variability across networks. This was done by calculating the maximum (S1/M1 and FPN) or minimum (DMN) single-trial amplitude of the BOLD response from each trial within a ± TR time window of the mean peak BOLD response. For each subject, we calculated the ratio of the single-trial amplitude standard deviation to the peak signal change measured from the mean BOLD response timecourse, and averaged over subjects. This metric expresses the proportional relationship between the trial-by-trial variability and the magnitude of a subject's average response.

CMRO 2 estimation
The following methodology was applied to the S1/M1 timecourses derived from the GLM analysis, and the S1/M1, DMN and FPN timecourses derived from ICA. The maximum signal change in the average BOLD and CBF responses were input to the Davis model (Eq. (1)) (Davis et al., 1998) to estimate the resultant change in CMRO 2 relative to baseline levels (CMRO 2 ) 0 for each subject: where the subscript 0 denotes the parameter value for the baseline period; α (Grubb coefficient) was chosen to be 0.2, in line with recent MR literature (Chen and Pike, 2010); and β, which reflects deoxyhaemoglobin concentration, was chosen to be 1.3 (Mark et al., 2011). M represents the maximum BOLD signal change, i.e. due to an increase in CBF, which causes complete oxygen saturation in venous vessels. M is dependent on field strength and TE (Chiarelli et al., 2007a) and is often calculated using a hypercapnic challenge (Gauthier et al., 2011). Since a hypercapnic challenge was not performed in this study, a range of M-values appropriate for grey matter was identified from the literature and normalised to the field strength (3 T) and TE (33 ms) used in this study, resulting in M values ranging from 6% (Gauthier et al., 2011) to 14.9% (Kastrup et al., 2002). The percentage change in CMRO 2 (Δ%CMRO 2 ) was calculated for each of these M values as well as for the mean value of M (10.45%). The coupling ratio (n) of Δ%CMRO 2 / Δ%CBF was calculated for each of these M values to ascertain the effect of variation in this parameter.

Results
In contralateral S1/M1, a significant increase in BOLD signal associated with a spatially coincident increase in CBF was found in response to MNS in all subjects using both ICA and GLM analysis. However, only ICA identified coherent spatial patterns of BOLD and CBF responses in the DMN and FPN, in addition to S1/M1. We first describe the spatiotemporal properties of BOLD and CBF responses in the S1/M1 network, and compare the Δ%CMRO 2 /Δ%CBF coupling calculated using the GLM and ICA methods of response localisation, before focusing upon the BOLD and CBF responses and Δ%CMRO 2 /Δ%CBF coupling in the DMN and FPN ICNs.

S1/M1: spatio-temporal correlation between BOLD and CBF response
Conventional GLM analysis using a boxcar regressor of consistent amplitude stimulation identified significant positive BOLD and CBF responses to MNS in the contralateral primary S1/M1 (Fig. 1A). Similar spatial patterns of both BOLD and CBF responses were identified in contralateral S1/M1 using ICA (Fig. 1B). A high degree of spatial similarity was observed between the group-level regions of positive BOLD and CBF response in contralateral S1/M1 (Figs. 1A & B, green voxels) using both GLM and ICA methods. A total of 1876 significant voxels were common between the BOLD-CBF conjunctions that were separately defined by GLM and ICA methods (Fig. 1C,  The temporal correlation between the entire BOLD and CBF timecourses was significant (p b 0.05) in 14/15 subjects for both GLM (R range = 0.08-0.67, median = 0.42) and ICA (R range = 0.04-0.45, median = 0.36) methods (note that the non-significant subject was different between GLM and ICA analysis). The temporal profile of the primary phase of the response was comparable between GLM and ICA methods, showing an increase in BOLD and CBF signals in all subjects, peaking 7-8 s post-stimulus onset (Figs. 1A & B). Group average peak response amplitudes were not significantly different between analysis methods for either BOLD (GLM analysis 0.54 ± 0.16%, ICA 0.47 ± 0.15%, p = 0.67), or CBF (GLM analysis 25 ± 10.1%, ICA 23 ± 9.3%, p = 0.81).

S1/M1: CMRO 2 estimation
The maximum signal change in the mean BOLD and CBF response timecourses extracted from the GLM and ICA data are plotted for each subject in Figs. 2A and B respectively. A significant positive correlation between maximum Δ%BOLD and Δ%CBF response amplitudes was found using the GLM analysis (R = 0.89, p b 0.001) and ICA (R = 0.82, p b 0.001). Data points were also observed to cross the CMRO 2 isocontours ( Figs. 2A & B) for both GLM and ICA data, indicating an increase in metabolism due to MNS.
The linear fit between Δ%CBF and Δ%CMRO 2 for GLM (n = 0.66 ± 0.02) and ICA (n = 0.67 ± 0.03) data are displayed in Figs. 2C and D respectively (M = 10.45). The effect of varying M upon the Δ%CMRO 2 / Δ%CBF coupling ratio (n) is shown in Fig. 2E. The values of n were found to be highly comparable between the GLM and ICA methods. GLM and ICA values of n were within error margins of each other across the range of M values investigated.
These results demonstrate that ICA can localize and extract both BOLD and CBF responses and allows the calculation of changes in oxygen metabolism to MNS which are equivalent to GLM methods. We further used timecourses extracted from ICA to investigate the spatial and temporal similarity of BOLD and CBF responses in the DMN and FPN (not identified by the GLM) and compared the Δ%CMRO 2 /Δ%CBF coupling ratios in these networks. Group Z-statistic maps and fMRI response timecourses of A) significant contralateral S1/M1 responses to MNS defined using a fixed-effects GLM; B) the contralateral S1/M1 component defined using temporal concatenation ICA, from BOLD (red) and CBF (blue) data. The conjunction of significant BOLD and CBF voxels is superimposed in green. Positive BOLD and CBF responses to MNS were observed in both GLM-and ICA-derived timecourses. The MNS period is indicated by the grey rectangle. Error bars represent the standard error of the mean over subjects. A high spatial overlap was observed between the BOLD-CBF conjunctions separately defined by GLM and ICA methods C).
high degree of spatial overlap was observed between BOLD and CBF ICNs for both the DMN (Fig. 3A) and FPN (Fig. 3B) in both of the group maps obtained using temporal concatenation ICA. Expressing this as the percentage of overlapping active voxels [=overlap / total active BOLD&CBF voxels, expressed as a percentage] resulted in the following values: DMN = 67%; FPN = 46% and single subject maps obtained using dual regression (DMN = 61 ± 14%, FPN = 53 ± 10%, mean ± std). fMRI response timecourses were obtained from each of the subjectspecific ICNs. The group average DMN timecourses showed a decrease in both BOLD and CBF signals below baseline levels in response to MNS (Fig. 3C), with the peak amplitude occurring between 10 and 12 s poststimulus onset. The group average FPN timecourse showed an increase in both BOLD and CBF signals above baseline levels (Fig. 3D), with the response amplitude reaching a peak between 7 and 9 s and being maintained for approximately 10 s. The temporal correlation of the entire BOLD and CBF signal timecourses was significant (p b 0.05) in 12/15 subjects for both DMN (R range = 0.03-0.57, median = 0.34) and FPN (R range = 0.01-0.51, median = 0.31) (note that the non-significant subjects were different between DMN and FPN).
We calculated a metric from the ICA timecourses to enable comparison of the single-trial response variability between networks. The group average ratio of the single-trial amplitude standard deviation to the mean BOLD response amplitude showed that the response variability was greater in the DMN (4.2 ± 3.2) and FPN (4.1 ± 2.8) than in S1/M1 (1.1 ± 0.8). This difference in the variability ratio was statistically significant (two-tailed, student's t-test) between the DMN and S1 (p = 0.003) and between the FPN and S1 (p = 0.001), but not between the DMN and the FPN (p = 0.92).

DMN and FPN: CMRO 2 estimation
The maximum signal changes in the mean BOLD and CBF responses that were extracted from ICA timecourses for each subject, are plotted against one another for the DMN and FPN in Figs. 4A and B, respectively. For the DMN, 12 subjects displayed negative Δ%BOLD and Δ%CBF mean response amplitudes. The mean amplitude of FPN Δ%BOLD and Δ%CBF responses was positive in 13 subjects. A significant positive correlation between maximum Δ%BOLD and Δ%CBF response amplitudes was found for both the DMN (R = 0.54, p b 0.05) and the FPN (R = 0.78, p b 0.01). Data points were observed to cross the CMRO 2 isocontours for both ICNs (Figs. 4A & B), suggesting that a decrease in metabolism occurred in the DMN and an increase in metabolism occurred in the FPN in response to MNS. The linear fits between Δ%CBF and Δ%CMRO 2 for the DMN (n = 0.8 ± 0.07) and the FPN (n = 0.65 ± 0.03) are displayed in Figs. 4C and D respectively (M = 10.45). Increasing the value of M had the effect of increasing the coupling ratio as shown in Fig. 4E. We observed that the values of n for the FPN were highly comparable to those measured in S1/M1 for all M values. However, higher n values were observed in the DMN than in the FPN (Fig. 4E), as evidenced by the non-overlapping error bars representing the standard deviation in the linear fit.

Discussion
In this study we used simultaneous BOLD and CBF measurements to investigate the BOLD, CBF and CMRO 2 responses to MNS in S1/M1 and the DMN and FPN ICNs. We localised spatially coincident increases in the BOLD and CBF signals in contralateral S1/M1, the cortical area that is directly driven by the afferent MNS input, using both ICA and conventional GLM analysis. BOLD and CBF timecourses derived from ICA were highly comparable to those extracted from a GLM-derived ROI in the preprocessed data. Consequently, we obtain very similar estimates of the change in CMRO 2 and the Δ%CMRO 2 /Δ%CBF coupling ratio (n) from using the two methodologies.
We subsequently used ICA to extract BOLD and CBF signals from the DMN and FPN ICNs, which were not identified using the standard GLM approach, presumably because of the non-canonical nature of their haemodynamic responses (Fig. 3). We found that the spatio-temporal patterns of signal fluctuations in the DMN and FPN are highly comparable between BOLD and CBF data. Using these responses, we estimated CMRO 2 changes in these networks during MNS, which suggested a larger metabolism/CBF coupling ratio in the DMN than was found in the FPN or in S1/M1. GLM vs ICA identification of S1/M1 responses ICA has become a popular technique for analysing fMRI data from both resting and task-based experiments. ICA provides a reliable, datadriven approach to decompose an fMRI data set into a set of maximally, spatially-independent maps and associated time courses (Beckmann and Smith, 2004;Calhoun et al., 2001b;McKeown et al., 1998). As the decomposition assumes temporally coherent signals, the neurophysiological, non-artefactual components can be said to represent networks of brain regions sharing common activity. ICA provides the opportunity to explore changes in fMRI signals which are either unexpected or undetected when using conventional GLM analyses that are constrained by an a-priori model relying on assumptions about the timing, shape and consistent amplitude of brain responses (Bartels and Zeki, 2004;Beldzik et al., 2013;Calhoun et al., 2001a,b;Eichele et al., 2008;Malinen et al., 2007;Schmithorst, 2005).
In agreement with previous work (Bagshaw and Warbrick, 2007;Calhoun et al., 2001a), we find a strong similarity between regional BOLD response timecourses extracted using ICA and those extracted from the GLM. We extend this work by showing good agreement between ICA and GLM methods for localising CBF sensorimotor cortex responses to MNS. Furthermore, estimates of the changes in S1/M1 oxygen metabolism induced by MNS were highly comparable between ICA and GLM analysis.
Previous studies using calibrated BOLD and CBF fMRI suggest that the Δ%CMRO 2 /Δ%CBF coupling ratio in sensorimotor cortex lies in the range n = 0.3-0.49 (Chiarelli et al., 2007a;Kastrup et al., 2002;Stefanovic et al., 2004). However it must be considered that these previous studies used either isometric contraction or finger tapping motor tasks, under the volitional control of the subject, which are not fully  FPN (B) spatial maps defined from BOLD (red) and CBF (blue) data show a high degree of similarity. The conjunction of significant BOLD and CBF voxels is superimposed in green. Group DMN and FPN spatial maps were identified using temporal concatenation ICA (G) and dual-regression was subsequently used to generate individual maps shown for five representative subjects (S1-S5). Group mean DMN (C) and FPN (D) responses to MNS for BOLD (red) and CBF (blue) data. The MNS period is indicated by the grey rectangle. Error bars represent standard error of the mean over subjects.
comparable to the passive MNS delivered here. In the current study we use fMRI responses to MNS in S1/M1 to calculate n = 0.66 ± 0.02 or n = 0.67 ± 0.03 for the GLM and ICA data, respectively (M = 10.45%). These values of n appear large compared with previous reports. However, the difference in n can be attributed to the difference in α and β values employed in the current study, which are in line with most recent literature (Chen and Pike, 2010;Mark et al., 2011). If the "classical" values (α = 0.38 and β = 1.3) are used, along with an M value of 9.5% taken from (Stefanovic et al., 2004) (adjusted for our field-strength and TE), then the resultant values, n = 0.49 ± 0.03 (GLM) or n = 0.50 ± 0.04 (ICA), are in agreement with previous reports for the motor cortex (Stefanovic et al., 2004). The close agreement of n values obtained for S1/M1 between ICA and GLM methodologies provides a validation of the use of ICA for estimating changes in CMRO 2 in response to stimulation. In general we find that there is no inherent advantage of using model-free ICA over conventional GLM methods for localizing and extracting responses from primary sensory cortical areas where the timecourse shape can be accurately predicted. Such responses commonly exhibit consistent amplitudes, canonical temporal profile and possess sufficiently high SNR to be extracted by either technique. However, recent work suggests that even simple stimulus paradigms recruit widespread brain regions (Gonzalez-Castillo et al., 2012; Mayhew et al., 2013), modulating activity in ICNs spatially and functionally distinct from the directly stimulated sensory cortex. These responses can exhibit complex response morphology, and therefore go undetected by a GLM analysis using a simple boxcar model of the underlying neuronal activity. In this scenario, ICA approaches become especially valuable for producing a more complete picture of whole-brain activity.
The influence of subtle changes in the activity of the attention, default-mode or saliency ICNs upon task processing has been demonstrated in a range of cognitive tasks Prado and Weissman, 2011;Sadaghiani et al., 2010;Spreng et al., 2010) indicative of an important role in supporting brain function. Therefore, enabling more detailed study of the physiological and metabolic changes that accompany this activity can help improve the understanding of these fundamental brain processes.

Task-positive and task-negative network responses
Group ICA was used to identify coherent patterns of both BOLD and CBF signal fluctuations in the FPN and DMN during the MNS paradigm that were not detected by the GLM. This result can be explained by a combination of the larger trial-by-trial variability we observed in the DMN and FPN response amplitude compared to that in S1/M1; the simple GLM boxcar model of stimulus timings; and deviations of the response timecourse from the morphology of the canonical HR used in the GLM analyses. The FPN fMRI responses display longer duration signal increases (N 20 s post stimulus cessation) whilst the time to peak (10-12 s) of the negative DMN responses are later than those typically modelled by the canonical HR. The deviations of the DMN and FPN from canonical responses illustrate the utility of ICA in identifying unpredictable brain activity. The amount of variability in the DMN and FPN responses is evident from the error bars on the group mean timecourses (Fig. 3). However it is clear that on average the FPN displays fMRI signal increases and the DMN displays fMRI signal decreases in response to MNS, although the average signal changes are smaller than those observed in S1/M1. These responses are consistent with the widely held description of the FPN and DMN as task-positive and task-negative networks respectively (Dosenbach et al., 2007;Gusnard et al., 2001;Hutchinson et al., 1999;Raichle et al., 2001), although it must be noted that this simplified conceptualisation does not hold in all experimental contexts (Spreng, 2012).

CMRO 2 /CBF coupling in the DMN and FPN
This study advances understanding of the metabolism-flow coupling in two of the brain's fundamental networks, the DMN and FPN. Our results strongly suggest that, for the MNS used here, the Δ%CMRO 2 /Δ%CBF coupling ratio is larger in the DMN than in the FPN and S1/M1. These findings are in contrast to a recent study which used an arithmetic calculation task to induce negative BOLD and CBF responses in DMN regions and reported that DMN n values were similar to those calculated in the occipital cortex and inferior frontal gyrus (IFG) regions, which responded positively to the same task, for M values ranging between 5.7 and 25% (Lin et al., 2010). Lin et al's upper estimate of n for the PCC node of the DMN was 0.62 ± 0.02 based on using M = 25% (Uludağ et al., 2004). This is very similar to their corresponding upper estimate of n = 0.62 ± 0.01 in the positively responding occipital cortex and IFG. However, the discrepancy between this study and the findings described here could arise from a number of experimental differences: 1) Lin et al. used a 4 T scanner and the traditional value of α = 0.38 (here we use 3 T and α = 0.2); 2) BOLD measures were obtained from an ASL sequence with TE = 17 ms in Lin et al. (2010), compared to our use of a DABS sequence optimised to measure concurrent BOLD (TE = 33 ms) and CBF (TE = 13 ms) signals; 3) the presence of higher physiological noise in non-background suppressed ASL data; and 4) Lin et al. also used a very different cognitively engaging task to induce strong DMN deactivations, compared to our passive, sensory MNS task.
The difference in the n values that we observe in the FPN and DMN suggests a difference in the Δ%CMRO 2 /Δ%CBF coupling that underlies the responses in these brain regions. However, when interpreting differences in Δ%CMRO 2 /Δ%CBF coupling and comparing these results to previous work in sensory cortex, it is important to consider the sensitivity of the Davis model to variations in the parameters α, β and M. As found with the Δ%CMRO 2 /Δ%CBF coupling in S1/M1, the calculated values of n in the DMN and FPN are higher than previously reported. However, it is difficult to make comparisons between our DMN and FPN measures and previous studies as, to the authors' knowledge, only Lin et al. have calculated n outside of primary sensory cortex (Lin et al., 2010). The use of recently updated parameter values in the current study can again partially account for this difference. Using the "classical" values of α = 0.38, β = 1.2 and M = 10.45% gave n = 0.46 ± 0.04/0.65 ± 0.06 for the FPN/DMN, which is similar to/larger than previous estimates in sensory cortex (Ances et al., 2008).
A difference between the coupling in the FPN and DMN was observed when using a range of α and β values (α = 0.15-0.45 and β = 0.9-1.4), suggesting that the values of α and β employed in this study are not the cause of the difference in n which we observe between the ICNs. We further tested the possibility that increased sensitivity to α and β in the task-negative DMN compared with the task-positive FPN could result in an over-estimation of CMRO 2 , and consequently n in the DMN. Equivalent n values for the two ICNs were only obtained (n = 0.64) when using α = 0.2 and β = 1.3 for the FPN and the combination of α = 0.3-0.42, β = 0.95-1.2 for the DMN. However no previous work supports the use of such divergent parameter values across brain regions.
The value of M employed in the Davis model also affects calculations of n (Chiarelli et al., 2007b). As hypercapnic calibration was not used to calculate subject specific M-values in this study, it is not possible to characterise the coupling ratios in the FPN and DMN definitively. Therefore, in this study we investigated the effect on the coupling ratio of using a range of M values derived from the literature, scaled to our field strength and TE. Fig. 4E strongly suggests that if the same M value can be used for both of these ICNs, then n is significantly larger for the FPN compared to the DMN, regardless of the exact value of M employed in the Davis model. We also note that even if a physiologically unrealistic value M = 100% is used, n is still found to be higher in the DMN than in the FPN.
Since the PCC and the mPFC nodes of the DMN are deeper brain structures with different composition and vasculature to that found in the upper cortical IPS and lPFC regions of the FPN, it is possible that the maximum BOLD response, and hence M, could vary between the two ICNs. It is further possible that M could vary between the different nodes of a single ICN that is comprised of widespread brain regions. However, the literature suggests the variability of M within grey matter is small compared to that between subjects (Gauthier et al., 2011). Although the value of M for the DMN and FPN has not been specifically calculated in previous studies, parameters characterising cerebral physiology across the whole brain have been estimated in several recent calibrated fMRI studies (Bulte et al., 2012;Gauthier and Hoge, 2012;Wise et al., 2013). These studies showed regional variability in M over the cortex and the data is suggestive of a higher value of M in the primary DMN node of the PCC than other grey matter brain regions (Gauthier and Hoge, 2012;Wise et al., 2013). The oxygen extraction fraction (OEF) is thought to be consistent between brain regions (Wise et al., 2013), therefore the value of M is closely influenced by the cerebral blood volume (CBV). As the PCC is known to exhibit relatively high CBV (Ito et al., 2003), a larger value of M could be expected for the PCC. However, Fig. 4E clearly shows that increasing M has the effect of increasing the value of n. Using a larger M for the DMN compared to the FPN would consequently cause a greater divergence of n across the two ICNs than using the same value of M.
To investigate the variability in coupling ratio between networks further, we analysed individual nodes of the DMN and FPN using the methods described above. For both BOLD and CBF group ICA data, the DMN component was manually divided into three ROIs (PCC, bilateral IPL and mPFC) and the FPN into bilateral IPS and bilateral lPFC ROIs. Dual regression was then used to extract BOLD and CBF response timecourses and Δ%CMRO 2 /Δ%CBF was calculated for each ROI. We found that n was largest in the PCC and IPL, however there was no significant difference in n between nodes in either network for the same value of M (Fig. S1 in Supplementary material).
In summary, we have not performed a hypercapnic/hyperoxic calibration, which might enable a more definitive value of n within the two networks to be calculated. However, previous investigations of the variation of M across the cortex indicate that the difference in n that we observe between the FPN and DMN is consistent with a difference in the CBF-metabolism coupling that underlies the responses in these brain regions, and that the coupling we report could be considered a lower limit given the effect of likely spatial variations in M.
The DMN has recently been shown to exhibit levels of resting-state aerobic glycolysis (linear regression of glucose utilization onto oxygen consumption) that are significantly higher than the brain average (Vaishnavi et al., 2010). Further work is required to elucidate whether this difference is also present during task-induced changes in DMN activity, but it may suggest that the DMN activity is supported by energy supply mechanisms different to those of other cortical areas, which is in line with the findings presented here.

Origin of altered metabolism-blood flow coupling in the DMN
Previous studies have demonstrated that n is not an absolute, fixed quantity even within a cortical area and can vary with experimental manipulations such as caffeine administration , stimulus duration (Lin et al., 2009) and attention (Moradi et al., 2012). The effect of experimental context may be especially important to the modulation of activity in the FPN and DMN, as the extent to which these networks are recruited by a task will vary depending on the particular experimental paradigm and its associated cognitive demands. For instance, the magnitude of DMN deactivation has been previously linked with the level of engagement in task performance (McKiernan et al., 2005;Pallesen et al., 2009;Singh and Fawcett, 2008) and FPN recruitment varies with attentional load (Adler et al., 2001;Lawrence et al., 2003;Nebel et al., 2005). Furthermore, the occurrence of functional interactions between the DMN and the FPN has been shown to predict and support the performance of certain tasks (Fornito et al., 2012;Prado and Weissman, 2011). It remains to be established what consequences such co-operative signaling and related activity between the two ICNs have upon Δ%CMRO 2 /Δ%CBF coupling ratios.
A difference in n between grey matter in visual cortex and subcortical thalamic nuclei has previously been reported (Ances et al., 2008) despite similar reported M values in both regions. Therefore, our higher n value for the DMN could be explained by a larger proportion of the DMN comprising deep cortical grey matter compared with the FPN. The coupling differences we observe between the two ICNs suggest a reduced responsiveness of CBF (compared with CMRO 2 ) in the DMN compared with the FPN, which may be due to differences in vascular biomechanics or control by the sympathetic nervous system (Ances et al., 2008). Additionally, a recent study has shown that CBF in the PCC node of the DMN was larger than the other brain regions at rest (Pfefferbaum et al., 2011). Furthermore, although CBF levels in this region decreased in response to a working-memory task, they remained higher than in the majority of other cortical areas investigated (Pfefferbaum et al., 2011), suggesting that the perfusion characteristics or neurovascular coupling of this region are perhaps not closely comparable to the rest of the brain.
Alternatively, an altered BOLD mechanism may underlie CMRO 2 variations in the task-negative DMN. The observation of similar coupling ratios between S1/M1 and the FPN leads to the possibility that the origin of the difference in n between the FPN and DMN is related to differences in the underlying neurophysiological basis of the positive fMRI response observed in the FPN compared with the negative fMRI responses observed in the DMN. Recently we have shown, using these data, higher Δ%CMRO 2 /Δ%CBF coupling in task-negative, ipsilateral S1/M1 compared with the task-positive contralateral S1/M1 (Mullinger et al., 2014). Negative BOLD responses to stimulation have been widely reported in sensory cortex as well as DMN areas and are thought to, at least partly, represent decreases in neural activity and metabolism (Kastrup et al., 2008;Schafer et al., 2012;Shmuel et al., 2006;Stefanovic et al., 2004). However, much still remains to be understood about the representation of inhibition in the BOLD signal, as recent work has demonstrated that increased inhibitory neuronal activity can result in a positive BOLD response (Enager et al., 2009;Pelled et al., 2009). Decreases in neuronal activity are known to arise through multiple mechanisms involving different contributions of decreased excitatory input and increased activity of inhibitory neuronal populations (Attwell et al., 2010;Boorman et al., 2010;Buzsáki et al., 2007;Cauli et al., 2004;Lauritzen et al., 2012;Logothetis, 2008;Moraschi et al., 2012;Schubert et al., 2008). However, the recruitment of different classes of neurons/neuronal-astrocyte interactions could result in different metabolic demands and different control of the vascular response  between the negative BOLD regions of the DMN compared with the positive BOLD regions of the FPN and S1/M1. This provides another plausible explanation for the difference in coupling ratios observed here and requires further investigation.

Conclusion
Our work demonstrates the importance of simultaneous BOLD and CBF measurements for relating coherent BOLD signal fluctuations in ICNs to physiological variables, such as the Δ%CMRO 2 /Δ%CBF coupling. Our results suggest that the default mode network may exhibit metabolism-flow coupling distinct from many other brain regions. Further investigation of the source of the observed metabolic differences between the DMN and FPN and S1/M1 cortical areas and extension into the study of other ICNs is important to improve understanding of the healthy brain's functional architecture and the significance of changes that occur in disease conditions.