Iterative analysis of cerebrovascular reactivity dynamic response by temporal decomposition

Abstract Objective To improve quantitative cerebrovascular reactivity (CVR) measurements and CO 2 arrival times, we present an iterative analysis capable of decomposing different temporal components of the dynamic carbon dioxide‐ Blood Oxygen‐Level Dependent (CO 2‐BOLD) relationship. Experimental Design Decomposition of the dynamic parameters included a redefinition of the voxel‐wise CO 2 arrival time, and a separation from the vascular response to a stepwise increase in CO 2 (Delay to signal Plateau – DTP) and a decrease in CO 2 (Delay to signal Baseline –DTB). Twenty‐five (normal) datasets, obtained from BOLD MRI combined with a standardized pseudo‐square wave CO 2 change, were co‐registered to generate reference atlases for the aforementioned dynamic processes to score the voxel‐by‐voxel deviation probability from normal range. This analysis is further illustrated in two subjects with unilateral carotid artery occlusion using these reference atlases. Principal Observations We have found that our redefined CO 2 arrival time resulted in the best data fit. Additionally, excluding both dynamic BOLD phases (DTP and DTB) resulted in a static CVR, that is maximal response, defined as CVR calculated only over a normocapnic and hypercapnic calibrated plateau. Conclusion Decomposition and novel iterative modeling of different temporal components of the dynamic CO 2‐BOLD relationship improves quantitative CVR measurements.

Reliably measuring quantitative CVR, however, is challenging due to various dynamic response components contained in the BOLD-CO 2 relationship. For instance, delay of the BOLD signal response to PaCO 2 (i.e., temporal delay) , varies for individual vessel response times and CO 2 arrival times. Contrarily, static CVR measurements -by using two different steady states, e.g. normocapnia und hypercapnia -only represents the maximal overall response (Vagal, Leach, Fernandez-Ulloa, & Zuccarello, 2009), thus leaving out important dynamics of the (patho-) physiological BOLD response.
To better model true CVR response, investigators (Blockley, Driver, Francis, Fisher, & Gowland, 2011;Donahue et al., 2015;Geranmayeh, Wise, Leech, & Murphy, 2015;Halani, Kwinta, Golestani, Khatamian, & Chen, 2015;Poublanc et al., 2013) have aimed to define the arterial arrival time of PaCO 2 by applying maximum Pearson productmoment correlation coefficient between CO 2 and BOLD signal response, assuming a constant rate of vascular dynamic response over the brain. In this instance, the delay inducing the maximum correlation is considered the arterial arrival time. The assumption of constant rate of vascular dynamic response over the brain has recently been challenged. Already for healthy vessels there is a different vascular response to CO 2 between subjects (Regan, Duffin, & Fisher, 2013;Regan, Fisher, & Duffin, 2014) as well as for white and gray matter vessels within a single subject (see Figure 4 from Bhogal et al., 2015). Besides, even at a similar rate of dynamic vascular response, the maximum correlation does not only incorporate arterial arrival time of CO 2 , it also integrates part of the vascular dynamic response function to PaCO 2 . Therefore, calculating CVR after maximum correlation may result in suboptimal separation of CO 2 arrival time and consequently a location-dependent erroneous CVR interpretation.
Using an iterative analysis, we model and separate different components of the dynamic CO 2 -BOLD relationship to acquire more detailed information about the dynamic vessel response times to CO 2 between different brain regions on a voxel-by-voxel basis. By redefining the CO 2 arterial arrival time, we are able to provide a more optimal dynamic CVR calculation for (patho-) physiological delayed arrival time for every voxel in the brain. Additionally, this method allows for BOLD MRI-based static CVR measurements.

| Data selection and acquisition
This study was approved by the cantonal ethics board of the Canton of Zurich (KEK-ZH-Nr. 2012-0427). Normal datasets from 25 subjects were selected from an ongoing prospective study on BOLD-MRI CVR at our institution. Exclusion criteria were the presence of severe cardiopulmonary disease, MRI contraindications, age <18 old or the inability or refusal to sign informed consent. For this analysis, two subjects with unilateral carotid artery occlusion

| Iso-oxic CO 2 stimulus
End-tidal partial pressure of oxygen (PetO 2 ) and carbon dioxide (PetCO 2 ) were precisely controlled and targeted using a preprogrammed sequence, executed by a custom built computer controlled gas blender (Prisman et al., 2008;RespirAct TM , Thornhill Research Institute, Toronto, ON, Canada) with prospective gas targeting algorithms described by Slessarev et al. (2007). Specific details and novelties of this technique as compared to other vasoactive stimuli have been described more specifically in a previous publication (Fierstra et al., 2013). The preprogrammed sequence consisted of an initial 100 s baseline P et CO 2 at 40 mmHg, after which a P et CO 2 step was increased to 10 mmHg above baseline for 80 s and a return to baseline for 100 s, before free breathing was restored. Iso-oxia was maintained at 100 mmHg (Figure 1 -green line).

| Spatial preprocessing
Anatomical and functional images were preprocessed using Statistical Parametric Mapping 12 (SPM 12, Wellcome Trust Centre for Neuroimaging, Institute of Neurology, University College London, UK). First BOLD images were realigned to the mean BOLD image.
The high-resolution MPRAge T1-weighted image was coregistered to the mean BOLD image and probability maps for gray matter, white matter, cerebrospinal fluid, skull, skin and air were generated.
F I G U R E 1 Flow Chart DTP calculations. Flow charts of the iterative DTP determination on a voxel-wise basis. The temporal decomposition uses a maximum of 15 iterations. DTP, Delay to Plateau Left and right hemispheres were manually segmented. For construction of normalized volumes, functional and anatomical maps were normalized to the MNI template using a non-linear transformation into MNI space. After removal of all voxels, excluding the white and gray matter voxels to decrease the partial volume effects of neighboring CSF voxels, the functional images were spatially smoothed with a Gaussian kernel of 8 mm full width at half maximum.

| Temporal preprocessing
The MR data sets were further preprocessed with custom MATLAB R2013b routines (The MathWorks, Inc., Natick, Massachusetts, United States; http://www.mathworks.com/). A low band-pass filter with a filter cut-off frequency of 0.125 Hz was applied to the MR data following Duffin et al. (2015) (Figure 2). Thereafter, the BOLD time series were temporally smoothed by a 16 dynamics (6%) local regression using weighted linear least squares and a secondorder polynomial model with assignment of lower weight to outliers (robust Loess method). This procedure did not induce any time shift in the data. Last, we linearly detrended the BOLD time series.
Additionally, the PetCO 2 trace was resampled to match the TR of the BOLD data.

| CO 2 arrival time (Temporal Delay)
We used three different types of CO 2 arrival times (Temporal Delay) to match the CO 2 to the BOLD time course: 1. a global CO 2 shift for the whole-brain (Delay global ), the most commonly used method .
For the first two types (Delay global and Delay maxcorr ), we implemented a linear least-square best fit of these CO 2 time course to the BOLD signal for maximum cross-correlation (Pearson Correlation). The lower and upper ranges for the delay determination were −10 and 50 TR (80 s). We did allow for a small negative cross correlation. An upper limit of 50TR was chosen since a delay exceeding 50 TR(100 s) was considered physiologically improbable. The delay with the maximum Pearson correlation was used to calculate the CVR maps, see section 2.8.
F I G U R E 2 CO 2 time evolution vs. BOLD signal in two illustrative voxels. For each processing method, the shifted CO 2 timeline (green) is plotted against the BOLD signal time series of two random voxels (one in the white matter: red, and one in the gray matter: blue) together with the linear fitting of the BOLD vs. CO2 scatter plot and the respective CVR map. In figure 2.1, the CO 2 time series correspond to CO 2 time course without time shift. figure 2.2 shows a single global CO 2 time shift for the whole-brain (Delay global ) and figure 2.3 correlation CO 2 time shift (CO 2 -BOLD) on a voxel-wise basis after maximum correlation (Delay maxcorr ); figure 2.4 shows the respective duration of the DTP and the DTB for both voxels with the dynamic (CVR corr ) and static (CVR stat ) CVR maps. The fitting of figure 2.4 is based on CVR stat , scatter plot which shows clearly removal of the transition phases between the two static states. The shifts of the CO 2 timeline are: Delay global : red: 7 TRs, blue: 7 TRs; Delay maxcorr : red: 20 TRs, blue: 5 TR; Delay corr : red 9 TR, blue: 1 TR. CO 2 , Carbon dioxide; BOLD, Blood-oxygen-level-dependent; CVR, cerebrovascular reactivity; DTP, Delay to Plateau; DTB, Delay to Baseline

| Delay to plateau
To characterize the dynamic response we used the ten-to-ninety (10% and 90%) rise time duration. Conceptually such a method is used to determine the transition duration between two steady states for instance in electrophysiology (Floyd, 2014). We applied a similar concept to determine Delay to Plateau (DTP), defined as the duration between the 10% (i.e., start) and 90% (i.e., stop) of the dynamic BOLD signal response to a rapid (pseudo square wave) P et CO 2 increase. The substantial difference with standard rise time calculations is that the BOLD signal either increases or decreases to a P et CO 2 increase. The calculations, done on a voxel-wise basis, are explained in the flow chart in Figure 1. The maximal number of iterations was 15.

| Corrected delay
Corrected Delay (Delay corr ) values were calculated as followed: where Volume 10% is the index of the BOLD volume where the BOLD signal increases by 10% -at the beginning of the DTP. Volume start is the BOLD volume index with the computer driven start of the CO 2 step change minus 10 TR. In determining Volume start , we assumed a negative of 10 TR to correct for any delay between the display of CO 2 on the computer and the actual CO 2 inspiration, For each patient, we determined the fastest responding voxel and corrected the other voxels accordingly. This allows for optimal inter-subject analysis. Consequently, Delay corr is a voxel-wise corrected temporal delay in seconds between these two BOLD volumes indexes, corrected with the delay of the fasted responding voxel.

| Delay to baseline
Delay to Baseline (DTB) is defined as the time between the 90% to 10% signal return to baseline after the CO 2 step change. DTB is comparable to the commonly used "decay time". However, while decay time is usually seen as the signal drop between fixed percentages, DTB includes calculations for either increasing or decreasing signal responses, that is return to baseline.
For calculating DTB, the dynamic response from the end of the step change in CO 2 , a similar approach was used as for the DTP. We first determined the 10% Volume Index in a similar fashion as we did the 90% Volume Index for DTP. Thereafter, we determined the first Volume Index with BOLD signal intensity >90%. Similar to DTP, this process was also reiterated. The amount of indexes between the 90% and 10% Volume Indices, multiplied by two, is determined to be DTB.

| CVR calculations
After shifting the P et CO 2 time course, P et CO 2 was regressed against the BOLD signal time course using a linear least square fitting to calculate CVR, defined as: where %ΔBOLD is the percent of BOLD signal change and ΔCO2 in mmHg is the CO2difference. For comparison, CVR calculations were done with BOLD time series versus the global shifted CO 2 time series and a shifted CO 2 time series on a voxel-per-voxel basis (Delay maxcorr ) to create CVR global and CVR maxcorr (See Figure 2; Poublanc et al., 2013;Geranmayeh et al., 2015).
Both CVRs and Delays were color coded and presented as an overlay on the high resolution T1-weighted image, only to include voxels passing the 0.9 probability threshold in the combined gray and white matter probability map.

| Corrected CVR and static CVR
Corrected CVR (CVR corr ) is calculated in a similar fashion as CVR maxcorr with the difference that Delay corr is used as the temporal voxel-wise CO 2 shift.
Static CVR (CVR stat ) differs from CVR corr as the DTP and DTB periods were excluded from the calculations to determine maximal CVR response by only using plateau phases (here isocapnic and hypercapnic plateau phases; See Figure 2). , Sobczyk, Battisti-Charbonney,  and Poublanc, Crawley, Sobczyk, Montandon, et al. (2015), have previously defined the approach to generate CVR reference atlases, with increased sensitivity for impaired CVR. We have created our reference atlas in a similar fashion.

| Reference atlases and Z maps
For CVR corr , CVR stat , Delay maxcorr , Delay corr , DTP and DTB, the mean and standard deviation over the 25 datasets were calculated to establish the reference atlases.
After creation of each reference atlas with standard deviation, we calculated the Z scores on a voxel-wise basis as followed : where value Z is the score for the aforementioned variable (e.g., CVR corr CVR stat , etc.), X is the voxel value for the given variable, μ is the mean value for the given variable of our reference atlas for that voxel and σ is the standard deviation for the given variable in that voxel.
Abnormal voxels were defined as: | z scores | > 2σ, which implies a t-test value smaller than 0.05. Z-maps were color-coded and include just voxels with values outside the normal variation range.

| Decomposition of dynamic temporal components
For our iterative analysis we separated the dynamic temporal components and named them as follows: CO 2 arrival time -temporal delay-, Delay to Plateau, Delay to Baseline. Figure 2 illustrates the concepts of different CO 2 arrival times and DTP and DTB for a gray matter and a white matter voxel of an exemplary subject.
As becomes apparent, Delay global is too short for deep white matter voxels (red - Figure 2.2), and too long for gray matter and subcortical voxels (blue - Figure 2.2). Delay global calculations often match a gray matter maximum correlation more optimally. The implementation of a voxel-wise delay results in improved CVR contrast mainly in the white matter.CO 2 time series shifted with Delay maxcorr in Figure 2.3 partly corrects the aforementioned limitations inherent to Delay global .
As Delay global calculations often match a gray matter maximum cor- Removing DTP and DTB from the dynamical considerably improves the fit which we termed CVR stat . Illustrative functional maps of two subjects can be found in Figure 3a,b. Furthermore we have included illustrative axial slices of DTP, DTB, Delay maxcorr , Delay corr , CVR corr and CVR stat for both subjects in the supporting information.

| Reference atlases
For the formation of the reference atlases for DTP, DTB, Delay maxcorr , Delay corr , CVR corr and CVR stat , a total number of 25 subjects (age 34.2 ± 11.8) were included. Subject characteristics and extended slices of all six reference atlases are given in the supporting information (see Table S1 and Figure S1).
The average whole-brain DTP and DTB of our reference cohort are 25.1 ± 3.9 s, respectively 29.3 ± 4.5 s.
The overall Delay maxcorr over the whole-brain is 20.  Figure S2). This renders Delay corr , in combination with a small F I G U R E 3 Effect of arterial arrival time on CVR calculations. (a-b) shows a single axial slice of 2 CVR maps (CVR global and CVR corr ) and a Delay corr map of two representative subjects (age 24 and 26 respectively). CVR is color-coded between −0.2% and 0.3% to increase CVR contrast and highlight changes. The Delay corr is presented as an overlay on a T1-weighted image and color-coded between 0-80 s. Although subtle, it is obvious that in areas with prolonged Delay corr negative CVR findings with CVR global are corrected after CVR corr calculation. Moreover, noticeable is that CVR contrast does not increase in all regions of the brain after Delay corr implementation. In areas in the white matter without prolonged CO 2 arrival time, the Delay global presents a Delay value more closely matching the Delay maxcorr then the Delay corr . This results in an increased correlation between BOLD and CO 2 time series, enhancing CVR contrast. However, as is shown in Figure 2, this CO 2 arrival time also incorporates part of the dynamic response, which results in erroneous CVR calculations.
standard deviation very sensitive to detect abnormal voxels. Overall CVR values for CVR corr are generally smaller than CVR stat (0.26 ± 0.11 vs. 0.33 ± 0.14).
From the reference atlas, it becomes apparent that the DTP and DTB differences between individual subjects are fairly small (on average 25.3 ± 4.0 vs. 28.8 ± 4.5 s).

| Functional maps of patient datasets
For two subjects with unilateral ICAO an exemplary BOLD-derived functional axial slice of each parameter in combination with an axial slice of the reference atlas and abnormality analysis is given in In subject 2, in the right sided areas exhibiting impaired CVR, Delay maxcorr is highly influenced by the dynamics of DTP and DTB, while posteriorly to the impaired CVR an increased Delay corr is also present.

Both subjects experience abnormal voxel increases in their time
constant parameters at the borders for the impaired CVR area. Subject 2 also shows abnormal decrease within the affected area.
For subject 2, after removing the dynamic components, CVR stat shows a stronger increase in the ACA territory on the right. This is strongly related to the increased DTP found in this patient. For the z-score map, however, in the ACA only a small correction is observed, whereas the abnormal region in the MCA territory is enlarged.
Both patients show extensive and numerous abnormal z-CVRscore regions in the contralateral "unaffected" hemisphere illustrating that a unilateral vasculopathy may affect the entire brain.

| DISCUSSION
We present a novel CVR modeling by decomposing iteratively different dynamic components of the CO 2 -BOLD relationship thereby acquiring specific information about the varying vessel response times and CO 2 arrival times -that is temporal delay -on a voxel-by-voxel basis. We have characterized a new optimal variant of the CO 2 arterial arrival time (Delay corr ), and introduced novel functional parameters,

DTP and DTB. This method by itself further improves quantitative
CVR measurements, and patient specific voxel-wise parameter probability can be expressed in statistical z scores by matching it to a reference atlas.

| Delay to plateau and delay to baseline
DTP characterizes the duration of the BOLD signal evolution in response to a pseudo square wave of CO 2 (rapid stepwise CO 2 increase of 10 mmHg), while DTB represents the dynamic return duration to the calibrated baseline from the end of the CO 2 pseudo square wave (rapid stepwise CO 2 decrease of 10 mmHg). DTP and DTB are modifications of the widely adopted -"rise time" and "decay time" neuroimaging parameters.
Our findings of different dynamic vascular responses in gray matter and white matter are in good agreement with findings of other groups using CO 2 as the vasoactive stimulus (Andrade et al., 2006;Bhogal et al., 2015;Rostrup et al., 2000). The subtle duration differences, seen between the subjects presented in Figure 2, have also been described by Transcranial Doppler studies investigating the dynamic response of healthy subjects to CO 2 in the MCA territory (Regan et al., 2013(Regan et al., , 2014. These studies also found a fast initial response and a delayed secondary response as we observed in both voxels in Figure 2. We also found similarities between our DTP and BOLD-time-to-peak measurements done by Donahue et al. (2013) in healthy subjects.

| CO 2 arrival time
Earlier studies determined the arterial arrival time of CO 2 or Delay maxcorr with maximum correlation calculations and showed an improved explained variance (Andrade et al., 2006;Donahue et al., 2015;Geranmayeh et al., 2015). This erroneously created the issue of BOLD signal increase prior to CO 2 arrival, which is highly improbable physiologically (see Figure 2.3).
Based on the start of the dynamic vascular response we discarded the concept of Delay maxcorr and created the term Delay corr . Delay corr is an improved surrogate of cerebral arterial arrival time of CO 2 . Figure 2 clearly presents Delay corr as the better arterial arrival time definition.
The optimization of CO 2 arrival time by Delay corr can be observed as the parallel increase of BOLD signal and CO 2 time series, as well as the parallel decrease in BOLD and CO 2 after the step-function CO 2 decrease. This confirms Delay corr primarily reflecting the arterial arrival time without the vascular dynamic response. Large differences in Delay corr can be seen between gray matter and subcortical white matter versus the deep white matter. Gray matter and subcortical white matter are quickly perfused by larger arteries with sufficient anastomosis while in deep white matter the vasculature consists of smaller end-arteries (Bhogal et al., 2014). For healthy subjects, the use of Delay global results in erroneous findings of steal phenomenon that are only partially corrected if the correct arrival time is applied (Figure 3). Thomas et al. (2013) showed that an adequate CO 2 stimulus will result in negative CVR due to high signal CSF displacement by vasodilatory response. This could explain negative CVR periventricular after correction with Delay corr .
Furthermore, Delay corr also shows a better correspondence BAT are marginally shorter. Indeed, the reference/start point, with ASL, is the time and location of water labeling, which can be used as a reference for signal changes measured in different regions in the brain. With our method, the starting point is the PetCO 2 trace-thus there is a delay from the time the gas is inhaled to the time that the trace is measured by the system with the BOLD signal potentially reacting in between. By using a small negative delay and adjusting the delay to the fasted responding voxel, we corrected for all these factors.
This minor prolongation of Delay corr relative to DSC/ASL BAT is most likely related to a limitation in the calculations due to TR shift and the TR relative discretization noise. Static CVR calculations with BOLD has been previously introduced by Poublanc, Crawley, et al. (2015). Their steady state CVR is determined by convolving CO 2 with a parametrized hemodynamic response function to best fit the BOLD time series on a voxel-wise

| Dynamic versus static CVR calculations
base. This will result in CO 2 time series, very similar to the BOLD time series, making all calculations more or less linear and therefore static. However, problems may arise when interpreting this steady state CVR. A controlled stimulus is used because of its potential to precisely alter PetCO 2 by prospecting end tidal targeting with high correspondence to PaCO 2 (Prisman et al., 2008). Therefore, modeling PetCO 2 values during the dynamic range to best fit the BOLD time series will result in a fictitious CVR, not corroborated by true underlying PaCO 2 values.
Another technique for static or time independent CVR calculations is implementation of a controlled ramp protocol (Mutch et al., 2012;. Under controlled CO 2 changes it can be argued that a ramp protocol with a slowly progressive CO 2 increase can induce vascular dilatation throughout the brain with enough delay to level out DTP. Therefore, a ramp protocol can result in a better time independent -static -CVR calculation without loss F I G U R E 4 Two illustrative patients. Axial slices of 4 time constant maps and two CVR maps of two patients with unilateral internal carotid artery occlusion. The center images reflect the reference atlases for each parameter. The four time constant maps include DTP, DTB, Delay maxcorr, and Delay corr. The time constant maps are color-coded between 0 and 80 s. The CVR maps shown are the CVR corr and the CVR stat and are colorcoded between −0.6% and +0.6% BOLD signal change per mmHg CO 2 . The two outer images are z-score maps for abnormality assessment. Only voxels surpassing −2σ (Blue) or +2σ(Red) are shown in this image. CVR, cerebrovascular reactivity; DTP, delay to plateau; DTB, Delay to Baseline of power. However, this is only justified if the ramp protocol is within linear range of the BOLD-CO 2 sigmoidal curve . Outside this linear range, the BOLD versus CO 2 evolution needs to be considered to calculate CVR meaningfully.

| Limitations
The correct DTP and DTB determination needs the use of a stepwise CO 2 change. A slower CO 2 change would increase DTP and DTB and not allow a correct impulse response function/rate determination. In case of non-stepwise PetCO 2 trace, no correction can be made using, for example, the rise and decay constants of the PetCO 2 , as the higher frequencies would not be present in the stimulus.
For these calculations PaCO 2 needs to be stable. Suboptimal control of the physiological parameters will result in either upward or downward CO 2 drifts. Consequently, these calculations cannot be done using a non-standardized CO 2 stimulus.
An overall limitation for CO 2 arrival time calculations with fMRI is the temporal discretization on the TR. Our TR of two seconds could be shortened using Multiband acquisition (Setsompop et al., 2012).
Last, the younger cohort included in the reference atlas could have confounded our z-score calculations, as age related CVR changes have been described (De Vis et al. 2015, Bhogal et al. 2016).

| Caveats of the current study
We do not determine or calculate the underlying changes in blood flow, nor measure cerebral flow directly. The determined components (i.e., DTP and DTB) are transient response durations that includes all factors of the dynamic vessel response, including change in cerebral blood flow. We characterize the vascular BOLD signal.
For this paper we aimed to focus on the best analysis method and therefore chose to use a uniform baseline, as the analysis method works independent of a subjects CO 2 /O 2 baseline. We acknowledge that a this uniform baseline instead of a physiological baseline has limitations regarding the clinical interpretation of the BOLD response, but suits the method usability demonstration aimed here.
Due to the noisy nature of the BOLD signal, DTP and DTB determination needed spatial and temporal smoothing. Similarly to others BOLD fMRI studies, the signal to noise ratio (SNR) in the WM is lower, due to the inhomogeneity of the coil array sensitivity. This low SNR might disturb the response and delays estimations of BOLD signal responses below the noise level. We, therefore, used a spatial 8 × 8 × 8 mm FHWM Gaussian low pass filter and a temporal robust Loess smoothing method with a 32 s sliding -shifted with each TR -window, with which a local regression using weighted linear least squares weights and a 2nd degree polynomial model is performed, assigning lower weight to outliers in the regression. The 16 time points used in this analysis are sufficient to provide a good estimate of the local signal evolution -to determining fitting parameters, and noiseto determine outliers.
The smoothing is performed a posteriori and "centered", so no additional delay is induced through the smoothing. This smoothing method results in overestimating of the very short DTPs and DTBs.
However, it also provides the possibility for robust DTP and DTB

calculations.
On the other hand, for the very long DTP and DTB, the accuracy of finding the true physiological DTP and DTB for a calibrated baseline CO 2 depends mainly on the sufficient length of the CO 2 increase and return to baseline values (after the stepchange). We have opted for a stepwise CO 2 change of 80 s since preliminary experiments confirmed such a period would be an adequate duration for most vessels to fully react to the CO 2 challenge while maintaining high tolerability of subjects to hypercapnia. We extended the CO 2 baseline to 50 TR (100 s) to allow tissue with potential CO 2 storage to have a larger DTB than DTP (i.e., to react slower to a decreasing CO 2 than during an increase). However, in areas of the brain with less vascular den- Comparison to longer CO 2 pseudo square wave changes will determine the optimal stepwise CO 2 duration needed to calculate DTP, DTB and CVR stat .

| CONCLUSIONS
Iterative decomposition and novel determination of different components of the dynamic CO 2 -BOLD relationship improves sensitivity and reliability of quantitative CVR measurements.

ACKNOWLEDGMENTS
The authors like to thank the University Hospital Zurich MRI technologists with their help in data acquisition.

CONFLICTS OF INTEREST
JA Fisher is the chief scientist of Thornhill Research Inc., Toronto, Canada (TRI), a spin-off company from the University Health Network Toronto. The company developed the RespirAct TM , which is currently dispensed as a non-commercial research tool around the world. J.