Insights into cerebral tissue-specific response to respiratory challenges at 7t: evidence for combined blood flow and co 2 -mediated effects


Cerebrovascular tone is constantly being modulated to ensure adequate supply of oxygen and glucose to brain tissues ( Paulson et al., 1990 ; Iadecola and Nedergaard, 2007 ). Control arterioles that branch from larger feeding arteries respond to local changes in metabolism, pH, or arterial CO 2 partial pressure (P a CO 2 ), work to regulate cerebral blood flow (CBF), ensuring proper delivery of nutrients and effective waste removal. This sensitive process, whereby dynamic changes in vascular tone modulate regional CBF in response to vasodilatory stimuli is defined as cerebrovascular reactivity [CVR; ( Fisher et al., 2018 ; Liu et al., 2018 )]. In recent years, CVR mapping has emerged as a robust clinical tool to understand the vascular physiology underlying conditions like stroke ( Geranmayeh et al., 2015 ), moyamoya disease ( Conklin et al., 2010 ), arterial stenosis ( Mandell et al., 2008a ), and traumatic brain injury ( Len et al., 2011 ; Ellis et al., 2016 ; Champagne et al., 2020 ), as well as the hemodynamic processes associated with normal aging of cerebral tissues ( De Vis et al., 2015 ; Bhogal et al., 2016 ).

Moving beyond amplitude-based assessments of vascular properties (i. e., CVR) advanced post-processing analyzes have emerged to characterize temporal markers associated with the CVR response ( Duffin et al., 2015 ; Poublanc et al., 2015 ; Donahue et al., 2016 ; van Niftrik et al., 2017 ; Champagne et al., 2019 ). These markers can provide quantitative profiles that aid in the explanation of differences between vasodilatory response across cortical and sub-cortical tissues ( Donahue et al., 2014 ; Thomas et al., 2014 ; Champagne et al., 2019 ) as well as pathological tissues with impaired or temporally altered hemodynamics due to collateralization ( Donahue et al., 2016 ). The general consensus is that temporal characteristics of the vascular response reflect a combination of factors including arterial transit time, blood redistribution, and vascular response speed, which may be partially separated using a combination of hypercapnic (HC) respiratory challenges to stress the cerebral vasculature and hyperoxic (HO) respiratory challenges that have previously been implemented to act as endogenous contrast agents via O2-mediated changes in deoxyhemoglobin ( Blockley et al., 2013 ; Champagne et al., 2019 ).

Through examinations of cerebrovascular architecture, it is known that cortical and deep gray-matter (GM) tissues have high arteriolar density ( Moody et al., 1990 ), which are fed by larger vessels branching off major cerebral arteries (i. e., pial, lenticulostriate, and choroidal). During HC, arteriolar smooth-vessel mediated dilation decreases peripheral resistance in cerebral arteries leading to an increase in CBF through capillary beds in accordance with Poiseuille’s law ( McDonald, 1974 ; Edvinsson et al., 2002 ). The bulk of these arterial resistance changes are in reaction to rising extravascular P a CO 2 that, in turn, drives changes in blood pH and local vasodilatory mechanisms ( Yoon et al., 2000 ; Tripp et al., 2001 ; Chesler et al., 2008 ). During MR imaging of HC, fractional changes in BOLD signal arise from the relationship between regional increases in blood flow and subsequent decreases in deoxy-hemoglobin concentration ([d-Hb]), which lengthens T 2 relaxation ( Ogawa et al., 1990 , 1993 ). Assuming this coupling between CBF and [d-Hb] remains constant, differences in local CO 2 accumulation, sensitivity and cross-membrane diffusion rate ( Thomas et al., 2014 ), as well as the possible redistribution (or “ vascular-steal”) effects between vascular networks ( Conklin et al., 2010 ; Fierstra et al., 2010 ; Bhogal et al., 2015 ), may all bias CVR measurements derived using standard methods, emphasizing the need to include additional temporal marker for the cerebrovascular response.

While HC modulates the BOLD signal via changes in perfusion, the BOLD contrast can also be manipulated when combined with non-vasoactive HO breathing challenges ( Blockley et al., 2013 ). As the arterial partial pressure of oxygen (P a O 2 ) increases, complete hemoglobin saturation at the lungs drives the dissolution of abundant O 2 into the blood plasma. This O 2 -loaded blood traverses the cardio-vascular network until it reaches the capillary beds where it diffuses toward cerebral tissues to meet metabolic demand. In terms of the BOLD contrast mechanisms, the increase in O 2 supply results in a reduced unloading of O 2 from arterial hemoglobin and a concomitant decrease in venous de-oxyhemoglobin content; this lengthens the blood T 2 relaxation time, creating a similar BOLD signal response as what is observed during HC. This O 2 based mechanism has been proposed as a means through which arrival time ( Champagne et al., 2019 ) and cerebral blood volume (CBV) ( Lu et al., 2003 ) can be assessed and assumes limited vasoconstriction and minor reduction in venous CBV ( Bulte et al., 2007 ; Mark and Pike, 2012 ; Xu et al., 2012 ). Irrespective of potential vaso-constrictive effects, the bolus of highly saturated blood formed during HO breathing can act as a non-invasive endogenous contrast agent for estimation of bolus arrival time ( Blockley et al., 2013 ; Liu et al., 2017 ; Champagne et al., 2019 ) since only minor reductions in blood velocity would be expected. Consequently, HO may be used to calibrate HC-derived CVR delays, and untangle physiological factors that contribute to the (patho)physiological mechanisms associated with diseases causing vascular impairments.

In MR-based methods that probe CVR using a gas delivery system ( Slessarev et al., 2007 ; Mark et al., 2010 ; Fisher, 2016 ), HC is typically induced via inhalation of CO 2 administered using a block or ramp respiratory challenge [see Figure 4 in Liu et al. (2018)]. Block-based HC is a stimulus where the end-tidal CO 2 (P ET CO 2 ) is rapidly increased from baseline in a boxcar function (typically +5–10 mmHg), which is assumed to evoke a dynamic CVR response ( van der Zande et al., 2005 ). In contrast, the RAMP-HC stimulus involves a gradual increases in the P ET CO 2 which is assumed to allow for more progressive changes in local P a CO 2 ( Ringelstein et al., 1988 ; Claassen et al., 2006 ; Battisti-Charbonney et al., 2011 ; Bhogal et al., 2014b ). Although the step design is commonly used to characterize temporal CVR delays, no studies to date has yet to explore tissue-based differences in delays from the vascular response induced during a ramp design, and how those may differ from CVR delays derived from block designs. Based on the assumption that local CO 2 levels are increased progressively during the ramp, it may be hypothesized that the slower arteriolar reactivity in response to decreasing pH allows for less blood flow redistribution and steal effects between vascular territories and tissues, that in turn, lead to shorter CVR delays. Additionally, given that HO-based BOLD mapping can be leveraged to estimate bolus arrival time, it may be hypothesized that a combined boxcar HC and HO respiratory stimulus can provide additional insight with respect to possible redistribution effects within local vessels, in the event that the region-specific lag times remain relatively increased compared to neighboring vascular territories.

In this study, we integrate HC- and HO-based respiratory challenges to gather insight about factors driving tissues-based differences in temporal delays of the CVR response. High spatial-resolution BOLD datasets were acquired at 7T using both BLOCK- and RAMP-HC paradigms to study the potential differences in CVR delays between dynamic and progressive local changes in CO 2 . Finally, a novel protocol using a combined HC and HO block, referred as the “ BOOST” protocol, was used to study the effect of bolus arrival time on tissue-based CVR delays, in contrast to CO 2 sensitivity.

Materials and Methods


We retrospectively surveyed our seven Tesla subject database to identify subjects in which both RAMP- and BLOCK-HC designs were acquired using the same MR imaging readout. We were able to identify a total of nine subjects. In a subset of seven of these subjects, an additional paradigm was identified in which a HO–HC block (herein termed “ BLOCK BOOST HC+HO”) was also administered. Age and gender information for the main (herein termed SAMPLE. 1) and sub-group (herein termed SAMPLE. 2) are provided in Table 1 .

TABLE 1 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 1

Subject demographics and end-tidal measures.

Respiratory Paradigms

Respiratory challenges were delivered using a 3rd generation RespirAct TM (Thornhill Research Inc., Toronto, ON, Canada) system in combination with a rebreathing face mask. The mask was taped to ensure an air-tight seal using Tegaderm transparent dressing (3M, St. Paul, MN, United States). For both groups, the BLOCK-HC protocol consisted of a 90s increase in hypercapnia, preceded, and followed by a pre- and post-baseline period, respectively ( Figure 1A ). In the case of the combined BLOCK BOOST HC+HO experiment, a simultaneous increase in arterial O 2 was implemented ( Figure 1A ). To ensure sharp transitions in arterial O 2 , the level of the O 2 increase was limited such that the gas delivery system could evoke arterial changes within several breaths in line with the transitions possible when increasing arterial CO 2 ( Figure 1A ). The progressive protocol (i. e., RAMP-HC) consisted of a baseline period followed by a 60 s hypocapnic period and a subsequent 5-min period of progressively increasing hypercapnia, ending with a return to baseline ( Figure 1A ). Ramped CO 2 protocols are well suited to determine the precise vascular reserve capacity in a particular region. For example, impaired CVR will manifest as an early plateau in the BOLD-CVR response as vessels exhaust their dilatory reserve capacity. The progressive nature of the ramped stimulus means that the subject only experiences the high arterial CO 2 levels required to elucidate mild impairments toward the end of the paradigm. In contrast, protocols based on BLOCK stimuli are suited to reveal dynamic aspects of the BOLD-CVR response and are imprecise in terms of revealing the dose-response profile of the CO 2 stimulus. For this reason, CVR information can be obtained using a lower stimulus magnitude that spares the subject from prolonged high-level exposure that can potentially cause discomfort. With the exception of one subject, all three paradigms were acquired on the same day during the same scan session. The baseline and achieved end-tidal CO 2 and O 2 values across all paradigms used in our study are reported in Table 1 . Typical end-tidal CO 2 and O 2 traces along with corresponding GM and WM ROI responses are shown for examples of each paradigm in Supplementary Figure 1 .

FIGURE 1 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 2

Summary schematic of the respiratory paradigms and imaging outputs. Sample data showing the RespirAct TM (Toronto, ON, Canada) end-tidal trace(A) for CO 2 (P ET CO 2 , mmHg; orange line and purple dot) and O 2 (P ET O 2 , mmHg; blue line and red dot) during the hypercapnic (HC) and hyperoxic (HO) gas challenges used in this study: BLOCK HC only (top), BLOCK BOOST HC+HO (middle) and RAMP HC only (bottom). The Blood Oxygen Level Dependent (BOLD) signal acquired simultaneously(B) was denoised using a Bayesian wavelet approach on a voxel-by-voxel basis and then extracted from the gray matter(C) to compute the optimized BOLD regressor using the Rapid Interpolation at Progressive Time Delays (RIPTiDe) method. The final BOLD regressor was then cross correlated against the original denoised BOLD signal(D) at each voxel to compute voxelwise lag maps (seconds)(E). This was repeated for each breathing challenges shown in(A), for each subject.

Magnetic Resonance Imaging

This study was approved by the medical research ethics committee of University Medical Center Utrecht (UMCU) and written informed consent was obtained from all subjects. The experiments were performed according to the guidelines and regulations of the WMO (Wet Medisch Wetenschappelijk Onderzoek). Subjects were scanned on a Philips 7 Tesla MRI system using a dual channel transmit coil in combination with a 32-channel receive coil. Image based (3rd order) shimming was performed. Respiratory challenges were administered throughout a multi-slice single-shot GE-EPI BOLD scan (flip angle: 90°, TR/TE 3000/25 ms, EPI/SENSE factor 47/3, reconstructed resolution: 1. 5 mm × 1. 5 mm, slice thickness: 1. 6 mm, no slice gap, FOV: 217. 6 mm × 192 mm, slices: 43, respiratory paradigm dependent scan duration: 507s for RAMP-HC and approx. 330s for BLOCK-HC/BOOST HC+HO).

Computation of Vascular Reactivity and Temporal Lags

Post-processing of the BOLD data consisted of temporal re-alignment [FSL: MCFLIRT; ( Jenkinson et al., 2012 )], segmentation [FSL: FAST; ( Zhang et al., 2001 )] and voxel-wise temporal de-noising using a Bayesian wavelet based approach ( Figure 1B ; symlet four wavelet, two levels, hard co-efficient threshold with level dependent noise estimation). End-tidal CO 2 and O 2 traces were resampled to the TR of the BOLD acquisition and aligned to the BOLD data based on the maximum correlation between the breathing trace and the individual mean GM signal. BOLD data and respiratory traces were then interpolated to a temporal resolution of 3000/8 ms (i. e., 8× oversampled) to account for temporal delays between slice acquisitions. The interpolated P ET CO 2 traces were used as the initial probe to generate and optimized BOLD signal regressor for a correlation-based temporal lag analysis [RIPTiDe; ( Frederick et al., 2012 ; Tong and Frederick, 2012 ); Figure 1C ]. Based on voxel-wise lag estimates, individual lag maps were generated for the BLOCK-HC, RAMP-HC, and BLOCK BOOST HC+HO respiratory paradigms ( Figure 1D ). P ET CO 2 traces as well as lag-adjusted traces were then used to generate CVR (Δ%BOLD/mmHg) and lag-adjusted CVR maps (seconds; Figure 2 ). This was accomplished by mapping the slope parameter obtained through linear regression of the lag adjusted P ET CO 2 trace against the voxel-wise BOLD signal. Finally, BOLD data were spatially normalized to the 1-mm MNI152 atlas via affine and non-linear bspline transforms using elastix [version 5. 0, ( Klein et al., 2010 ; Shamonin et al., 2014 )]. Transformation matrices were applied to CVR and lag-maps (using transformix) from each design to facilitate group region of interest analysis in segmented GM and WM tissues. Voxels containing less than four data points were thresholded from the MNI averaged CVR and lag maps.

FIGURE 2 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 3

Corrected cerebrovascular reactivity mapping based on RIPTiDe lag times. Sample data showing BOLD EPI axial slices from the BLOCK-HC challenge in a single subject(A) used to compute the RIPTiDe lag maps(B). The slice number is shown in green. In(C), the matching voxelwise correlation coefficient (r) map (axial, 20) shows the fit between the BOLD signal and optimized RIPTiDe regressor. The original and lag-corrected cerebrovascular reactivity (CVR) maps are shown in(D) and(E), respectively, with a specific look at deeper gray- and white-matter structures, to reflect the improvement in CVR following the lag time correction. Note that cerebrospinal fluid was removed from all maps to improve visual representation of the data. BOLD, Blood Oxygen Level Dependent, HC, hypercapnia, RIPTiDE, Rapid Interpolation at Progressive Time Delays.

Comparisons of Vascular Reactivity and Temporal Delays Across Tissues and Paradigms

Segmented tissue-based probability maps were binarized from the 1-mm MNI152 atlas using a voxel-based threshold set to 50 and 90%, for the GM and WM, respectively. The tissue-specific masks were applied to the group-averaged CVR and lag maps to extract the parameter distribution for each respiratory design, along with mean and standard deviation (MATLAB, version 2019a, The Mathworks, MA, United States). This procedure was repeated for both SAMPLE. 1 and SAMPLE. 2, in order to assess the repeatability of the CVR and delay measurements from the BLOCK- and RAMP-HC. Tissue-based histograms for each average map were also assessed to look at the distribution of the CVR lags, between GM and WM.


End-Tidal Measurements

All end-tidal measurements for each respiratory design are reported in Table 1 . On average, participants were exposed to an 8 ± 1 mmHg increase in P ET CO 2 during the BLOCK-HC, while maintaining changes in P ET O 2 fluctuations well below physiological threshold. Similar increases in P ET CO 2 were obtained during the BLOCK-BOOST (SAMPLE. 2), with concurrent increases in P ET O 2 averaged to 159 ± 44 mmHg, across subjects. During the ramp, larger increases in P ET CO 2 were obtained, ranging between [20. 8–22] ± 6 mmHg, with minimal changes in P ET O 2 .

Cerebrovascular Reactivity and Temporal Delays

All parameters described below are summarized in Table 2 . Average CVR measurements for both GM and WM tissues were in agreement between the BLOCK- and RAMP-HC paradigms. In comparison, CVR measurements were higher during the BLOCK-BOOST, although relative differences between the GM and WM were maintained. Qualitatively, similar observations can also be made based on whole brain assessment of lag-corrected CVR ( Figure 3 ).

TABLE 2 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 4

Gray- and white-matter average parameters for each stimulus.

FIGURE 3 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 5

Average lag-corrected cerebrovascular reactivity maps for each respiratory challenge.(A) The anatomical axial reference slices for the MNI template are shown with the slice number in green. The matching sagittal frame showing the referenced frames is displayed on the left. The group-averaged lag-corrected cerebrovascular reactivity (CVR) maps are displayed for each respiratory challenge in(B–D), representing the BLOCK-HC, BLOCK BOOST HC+HO, and RAMP-HC, respectively. A schematic of each respiratory design is shown on the left with the targeted end-tidal traces for CO 2 (yellow) and O 2 (cyan blue), as a reference.

While GM delays were similar between the BLOCK- and RAMP-HC (within ±2 s), longer WM lag times were observed during the RAMP-HC, when compared to the BLOCK-HC design ( Table 2 and Figure 4 ). This was repeatable across both samples studied ( Figure 4 ). In comparison to both HC designs, shorter GM and WM delays were observed during the BLOCK-BOOST HC+HO protocol, although the relative ∼10 s offset in the GM versus WM lag times was maintained.

FIGURE 4 Insights Into Cerebral Tissue-Specific Response to Respiratory Challenges at 7T: Evidence for Combined Blood Flow and CO 2 -Mediated Effects Picture 6

Summary of the tissue-based distribution of lag parameters for each respiratory design.(A) The anatomical axial slices for the MNI template providing a reference for the images presented.(B–F) The group-averaged lag maps (seconds) are displayed for each respiratory challenges, sub-divided based on the sample used to compute the mean image (SAMPLE. 1, N = 9, top, B–C ; SAMPLE. 2, N = 7, bottom, D–F ). The cumulative percent frequency (normalized to 100%) for the distribution of lag (seconds) is shown for each tissue which was extracted using the gray- (black) and white- (gray) matter mask displayed in(B), bottom right corner. A dotted red line was added to each histogram(B–F) at 40 s, for reference and comparison.

The r 2 fitting coefficient was similar for the BLOCK-HC and -BOOST, within both tissues ( Table 2 ). Greater fitting estimates were observed in the RAMP-HC, compare the block designs, in both the GM and WM masks ( Table 2 ).


In this study, we conducted a comprehensive analysis of the response to respiratory gas challenges in GM and WM tissues at 7T, as a way to advance our understanding of the effects that drive temporal delays in the BOLD-CVR response of healthy tissues. Findings from this study were three-fold: (1) In comparison to WM delays during the BLOCK-HC, longer WM delays were observed during the RAMP-HC suggesting that the progressive nature of the ramp stimulus may predispose the WM vasculature to a slower CVR response. We propose that this is a result of the delayed increase in the local intravascular CO 2 gradient between blood and tissue compartments. This delayed gradient reduced the driving force facilitating the diffusion of CO 2 into the tissues, which in turn, delayed the vasodilatory response. (2) The addition of the HO step within the HC boxcar design for the BLOCK-BOOST HC+HO protocol induced a global decrease in GM and WM delays despite maintaining a ∼10 s offset between the tissues, when compared to BLOCK-HC. We postulate that the vasoactive effect of CO 2 within tissues remained constant and that the shortening of the delays during BOOST may have been driven by blood arrival effects (i. e., non-vasoactive) arising from the endogenous O 2 contrast agent. (3) Comparable magnitude for lag-corrected CVR measurements were reported between the STEP-HC and RAMP-HC protocol, emphasizing that the step design is an appropriate tool to assess CVR markers, despite previous literature suggesting that a ramp-like stimulus may better model the sigmoidal relationship between changes in P ET CO 2 and CBF, driving BOLD changes in signal ( Bhogal et al., 2014a ). Altogether, these findings support the hypothesis that differences in temporal components of CVR between vascular networks reflect the culminative effect of CO 2 sensitivity (and/or CO 2 diffusion rate) in local vessels, in addition to blood flow redistribution and steal effects, as previously described. Moreover, these results suggest that the addition of a BLOCK-BOOST HC+HO paradigm within clinical settings can provide insights into whether diseases causing changes in CVR do so by way of severe blood flow redistribution (which would increase blood arrival time and inflate BLOCK-BOOST HC+HO delays), or alterations in vascular properties within the vessels that could impair CO 2 diffusion across WM tissues (i. e., atherosclerosis, amyloid angiopathy, aging-related fibrosis or stiffening of vessels; this would significantly increase the marginal gap in relative delay estimates between BLOCK-HC and BLOCK-BOOST HC+HO).

The GM and WM CVR values reported in this study are in line with previous studies looking at HC designs ( Bhogal et al., 2014b , 2015 ) at higher spatial resolution. That said, the WM CVR values are higher than ones reported using 3T imaging ( Thomas et al., 2014 ; Sam et al., 2016 ), likely due to a combination of factors including the higher BOLD CNR from 7T imaging, as well as the incorporation of temporal delay correction within the CVR computation, which has been shown to improve CVR estimates in healthy controls ( Duffin et al., 2015 ; Poublanc et al., 2015 ; Donahue et al., 2016 ; van Niftrik et al., 2017 ; Champagne et al., 2019 ). In this study, lag-corrected WM CVR was found to be lower in the WM, across all breathing paradigms. This is consistent with existing literature ( Brian, 1998 ; Rostrup et al., 2000 ; Mandell et al., 2008b ; Thomas et al., 2014 ) suggesting that despite temporal corrections for delays in the vascular response, WM tissue may have a lower vasodilatory capacity, at least in part, due to the longer time for extravascular CO 2 levels to build within the local vasculature, in response to global increases in arterial CO 2 content. As originally proposed in Thomas et al. (2014)(see Figure 3 ), the delayed intravascular build-up in CO 2 may be a function of the lower CBF within the WM ( Duvernoy et al., 1981 ), such that, changes in arterial content of CO 2 brought to the WM tissues per unit of time are comparatively lower than those for GM. This in turn increases the timespan required within which the intravascular-tissue CO 2 gradient can drive trans-membrane diffusion of CO 2 toward tissues, and modulate the local increases in extravascular CO 2 (which mediates the majority of the HC-induced vasodilation).

This theory is further supported by our observation of longer delays in the WM during the RAMP-HC respiratory paradigm, in comparison to the BLOCK-HC. During a step HC paradigm, the large and rapid increase in arterial CO 2 sets up a strong CO 2 gradient which remains high throughout the duration of the boxcar stimulus. This effect rapidly saturates GM vessels (due to early arrival time of GM blood) and continues to flow through penetrating arteries toward downstream WM vasculature. This large increase in intravascular CO 2 drives the concentration gradient along which CO 2 will diffuse to induce the vasodilatory response. In other words, a high gradient leads to a faster and stronger dilatory response whereas a low gradient may induce temporal delays irrespective of blood arrival time; particularly in regions having low overall CBV such as the periventricular WM. The rapid and dynamic nature of the step also provides a setting for physiological steal mechanisms and redistribution of blood flow effects ( Sobczyk et al., 2014 ; Bhogal et al., 2014a ; Poublanc et al., 2015 ; Champagne et al., 2019 ), by which hyper-sensitive areas may respond quicker to rising levels of arterial CO 2 , and inflate temporal delays in regionally closely bound vascular territories. This may be especially significant for pathological tissues that show impaired or temporally altered hemodynamic properties due to vascular collateralization ( Donahue et al., 2016 ).

In comparison to the BLOCK-HC, however, the RAMP-HC is assumed to be a progressive stimulus during which the vascular response has sufficient time to equilibrate with changing levels of arteriol CO 2 content (and the corresponding blood-tissue CO 2 gradient). Previously this was thought to lead to a more balanced response across brain regions through minimization of the dynamic redistribution (or steal) effects. This mechanisms may explain why GM delays during the RAMP were slightly higher than during the boxcar design (even though this difference was small; within ±2 s), as the slowly rising ramp stimulus allows for the vasculature to respond progressively in healthy tissues, in contrast to the dynamic STEP design which induces a rapid response in the GM vessels. In other words, healthy GM tissues may cope appropriately to rapid changes in arterial CO 2 associated with the BLOCK-HC, resulting in relatively shorter response time (and thus, CVR delays). This warrants further exploration in clinical populations where blood flow redistribution effects are much more likely (i. e., Moyamoya disease ( Conklin et al., 2010 )).

Contrary to our original hypothesis, our results showed that WM CVR delays were longer during the RAMP-HC, in comparison to the BLOCK-HC. We postulate that this occurred as a result of slowly increasing intravascular CO 2 associated with the ramp stimulus, paired with limited blood flow to the WM. Together, these effects organically delayed the build-up of higher CO 2 gradients within WM tissues, leading to an increase in the time required to evoke the expected vascular response in the WM. This is further exaggerated by additional steal effect from GM regions that are more sensitive due to higher vascular density (and therefore more CO 2 diffusion even considering lower CO 2 gradients). The response in WM increases notably toward the end of the stimulus paradigm, once the ramp reaches higher P ET CO 2 targets ( Figure 1A ) that setup stronger CO 2 driving gradients. These findings support the original hypothesis put forward in Thomas et al. (2014) providing an explanation for the higher CVR delays in WM within healthy brain tissues, which again, may be exacerbated in settings where regional blood flow and/or CO 2 sensitivity is disturbed.

Building on the findings from the comparison of the BLOCK- and RAMP-HC, the addition of the HO step during the HC challenge, using the BLOCK-BOOST HC+HO protocol, resulted in an overall decrease in the magnitude of delays for both the GM and WM. Despite those changes, however, the relative relationship between GM and WM was maintained during the BOOST protocol, showing a ∼10 s offset in delays between the two cerebral tissues. This supports the findings described above in that despite correction for arrival time of the stimulus into local tissues (signaled via the non-vasodilatory contrast O 2 agent), vasodilatory response delays in the WM were consistently higher than in the GM, suggesting that this phenomenon may result from slower CO 2 gradient build-up needed to drive diffusion into tissues. In other words, by using the BOOST protocol to compute CVR delays in clinical settings, in comparison to the BLOCK-HC delays, a more in-depth interpretation temporal response differences is possible, which teases effects related to the arrival time of the stimulus (where the magnitude of delays for BLOCK-HC and BOOST would be similar) versus local impairments in CO 2 sensitivity within WM tissues (which would inflate the ∼10 s offset gap between the GM and WM characterized between the two methods). The vasoactive mechanisms of CO 2 versus O 2 may differ and it is conceivable that their interplay could influence lag estimates. The elevated CVR resulting from the BOOST HC+HO stimulus when compared to the BLOCK-HC suggests that either the vasodilatory effects of the CO 2 overwhelmed any constrictive influence of O 2 , or the effect of increased S v O 2 increased the BOLD signal beyond any HO mediated flow reductions. However, considering the relatively small O 2 change applied in our experiments, we don’t expect appreciable constriction. Finally, the HO-induced BOLD signal is CBV-weighted, rather than reflecting CVR. This means that taken alone, the BOOST HC+HO paradigm may overestimate CVR since the HO-induced signal contribution is not explicitly accounted for. Therefore, although the BOOST HC+HO paradigm may provide additional temporal information about possible pathophysiological mechanisms, a standard CO 2 challenge is still required for accurate CVR magnitude measurements. Disentangling respective components of the BOOST HC+HO response remains interesting for future work.

Beyond the quantitative comparison of delays across the GM and WM, results from this study show that the distribution of CVR values across respiratory paradigm were repeatable, even after adding subjects from SAMPLE. 1 to SAMPLE. 2. This emphasizes that voxel-based CVR measurements are consistent within the brain of healthy controls, allowing for the clinical use of such biomarker as a tool to provide reference-based voxelwise analysis of vascular impairments in patients ( Sobczyk et al., 2015 ). Similarly, the consistency of the CVR delay distribution across the samples, between respiratory designs and within the GM and WM, supports that temporal analyses of the CVR response may too be used as a clinical tool to understand pathophysiological mechanisms associated with vascular diseases of the brain ( Donahue et al., 2016 ; Juttukonda and Donahue, 2019 ).

Despite the novelty of the findings presented in this study, some limitations must be acknowledged and recognized as opportunities for future studies. First, the study employed a retroactive method to analyze datasets that were scanned as part of branching research studies. This limited our ability to synchronize and put together complete datasets for all subjects used in the analysis. Moving forward, a greater cross-sectional research design may be employed to further assess the use of the BLOCK-BOOST HC+HO protocol, and its variability across a larger sample size. Furthermore, only healthy subjects participated, limiting the external validity of these findings and preventing our ability to predict whether possible changes in arrival time, or CO 2 sensitivity, as a result of vascular impairments, may be significant enough to be detected using the proposed methods. Future work should therefore consider implementing these analyses in clinical population with known vascular diseases, in order to assess the potential utility of such combined approach. Specifically, clinical studies may consider studying the effect of pathological disease mechanisms on the local diffusion rate of CO 2 into tissues, which could be tested using a combination of the BLOCK-HC and BLOCK BOOST HC+HO. It may be hypothesized that local damages to the endothelium may manifest as increases in the discrepancy between the GM and WM CVR lag times observed under each respiratory challenge, while blood arrival remains relatively unchanged (and thus delays would still be shorter under BLOCK BOOST HC+HO). Finally, the magnitude of the delivered BLOCK-HC change in P ET CO 2 was smaller than the maximum peak change in P ET CO 2 reached during the RAMP-HC, which may impact the values extracted for delays across tissues, given that the greater stress on local vasculature could promote a micro-environment within which blood flow redistribution effects are more apparent across regions with differences in hemodynamic capacity. Thus, future work may also consider examining the effect of changing the magnitude of the step stimulus and exploring whether larger CO 2 gradient changes the voxelwise lag maps.


As highlighted in this study, tissue-based differences in CVR and temporal markers for the response to hypercapnia are governed by a number of factors that reflects a compounded effect dependent on arterial arrival time, CO 2 sensitivity and CO 2 diffusion rate, blood flow redistribution, and steal effects, and vascular response speed. Here, we provide evidence showing that differences in the temporal components of CVR may be influenced by intravascular CO 2 gradients in local vessels, which, in addition to membrane permeability, and blood flow effects, determines the rate of diffusion within local tissues that then drives the vasodilatory response. The proposed mechanisms presented in this study, in support of Thomas et al. (2014), may provide a partial explanation for the slower WM response to hypercapnia in comparison to GM, in healthy tissues. Finally, the results presented suggest that the implementation of a BLOCK BOOST HC+HO respiratory challenge may provide additional insight about possible pathophysiological mechanisms underlying vascular diseases based on whether they are driven by severe blood flow redistribution, or alterations in vascular properties within the vessels that would affect trans-membrane CO 2 diffusion into tissues. Moving forward, the combination of a BLOCK-HC and BLOCK-BOOST respiratory paradigms may be used to untangle these factors driving CVR and the speed of response, as a tool to help to improve the diagnosis, prognosis and management of patients with vascular brain diseases.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by WMO (Wet Medisch Wetenschappelijk Onderzoek). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

AB was responsible for all the data collection and data pre-processing. AC ran the data analysis and wrote the manuscript. AB and AC edited and reviewed the manuscript and prior to submission. Both authors contributed to the article and approved the submitted version.


This research was funded by an NWO talent scheme grant (VENI) awarded to AB.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


The authors would like to acknowledge collaborational support from H. Hoogduin, M. E. P Philipens and J. C. W. Siero from the University Medical Center Utrecht in the Netherlands, as well as N. S. Coverdale and D. J. Cook from Queen’s University in Kingston, Canada.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www. frontiersin. org/articles/10. 3389/fphys. 2021. 601369/full#supplementary-material


Battisti-Charbonney, A., Fisher, J., and Duffin, J. (2011). The cerebrovascular response to carbon dioxide in humans. J. Physiol. 589, 3039–3048. doi: 10. 1113/jphysiol. 2011. 206052

Bhogal, A. A., De Vis, J. B., Siero, J. C. W., Petersen, E. T., Luijten, P. R., and Hendrikse, J. (2016). The BOLD cerebrovascular reactivity response to progressive hypercapnia in young and elderly. Neuroimage 139, 94–102. doi: 10. 1016/j. neuroimage. 2016. 06. 010

Bhogal, A. A., Philippens, M. E. P., Siero, J. C. W., Fisher, J. A., Petersen, E. T., Luijten, P. R., et al. (2015). Examining the regional and cerebral depth-dependent BOLD cerebrovascular reactivity response at 7T. Neuroimage 114, 239–248. doi: 10. 1016/j. neuroimage. 2015. 04. 014

Bhogal, A. A., Siero, J. C. W., Fisher, J. A., Froeling, M., Luijten, P., Philippens, M., et al. (2014a). Investigating the non-linearity of the BOLD cerebrovascular reactivity response to targeted hypo/hypercapnia at 7T. Neuroimage 98, 296–305. doi: 10. 1016/j. neuroimage. 2014. 05. 006

Bhogal, A. A., Siero, J. C. W., Fisher, J. A., Froeling, M., Luijten, P., Philippens, M., et al. (2014b). NeuroImage investigating the non-linearity of the BOLD cerebrovascular reactivity response to targeted hypo / hypercapnia at 7 T. Neuroimage 98, 296–305. doi: 10. 1016/j. neuroimage. 2014. 05. 006

Blockley, N. P., Griffeth, V. E. M., Germuska, M. A., Bulte, D. P., and Buxton, R. B. (2013). An analysis of the use of hyperoxia for measuring venous cerebral blood volume: comparison of the existing method with a new analysis approach. Neuroimage 72, 33–40. doi: 10. 1016/j. neuroimage. 2013. 01. 039

Brian, J. E. (1998). Carbon dioxide and the cerebral circulation. Anesthesiology 8, 1365–1386. doi: 10. 1097/00000542-199805000-00029

Bulte, D. P., Chiarelli, P. A., Wise, R. G., and Jezzard, P. (2007). Cerebral perfusion response to hyperoxia. J. Cereb. Blood Flow Metab. 27, 69–75. doi: 10. 1038/sj. jcbfm. 9600319

Champagne, A. A., Bhogal, A. A., Coverdale, N. S., Mark, C. I., and Cook, D. J. (2019). A novel perspective to calibrate temporal delays in cerebrovascular reactivity using hypercapnic and hyperoxic respiratory challenges. Neuroimage 187, 154–165. doi: 10. 1016/j. neuroimage. 2017. 11. 044

Champagne, A. A., Coverdale, N. S., Fernandez-Ruiz, J., Mark, C. I., and Cook, D. J. (2020). Compromised resting cerebral metabolism after sport-related concussion: a calibrated MRI study. Brain Imaging Behav. doi: 10. 1007/s11682-019-00240-2 [Epub ahead of print].

Chesler, M., Brien, J. O., Zappala, A., Cicirata, F., and Barrio, L. C. (2008). Regulation and Modulation of pH in the Brain. Brain 83, 1183–1221. doi: 10. 1152/physrev. 00010. 2003

Claassen, J. A. H. R., Zhang, R., Fu, Q., Witkowski, S., and Levine, B. D. (2006). Transcranial Doppler estimation of cerebral blood flow and cerebrovascular conductance during modified rebreathing. J. Appl. Physiol. 102, 870–877. doi: 10. 1152/japplphysiol. 00906. 2006

Conklin, J., Fierstra, J., Crawley, A. P., Han, J. S., Poublanc, J., Mandell, D. M., et al. (2010). Impaired cerebrovascular reactivity with steal phenomenon is associated with increased diffusion in white matter of patients with moyamoya disease. Stroke 41, 1610–1617. doi: 10. 1161/STROKEAHA. 110. 579540

De Vis, J. B., Hendrikse, J., Bhogal, A., Adams, A., Kappelle, L. J., and Petersen, E. T. (2015). Age-related changes in brain hemodynamics; a calibrated MRI study. Hum. Brain Mapp. 36, 3973–3987. doi: 10. 1002/hbm. 22891

Donahue, M. J., Faraco, C. C., Strother, M. K., Chappell, M. A., Rane, S., and Dethrage, L. M. (2014). Bolus arrival time and cerebral blood flow responses to hypercarbia. J. Cereb. Blood Flow Metab. 34, 1243–1252. doi: 10. 1038/jcbfm. 2014. 81

Donahue, M. J., Strother, M. K., Lindsey, K. P., Hocke, L. M., Tong, Y., and deB Frederick, B. (2016). Time delay processing of hypercapnic fMRI allows quantitative parameterization of cerebrovascular reactivity and blood flow delays. J. Cereb. Blood Flow Metab. 36, 1767–1779. doi: 10. 1177/0271678X15608643

Duffin, J., Sobczyk, O., Crawley, A. P., Poublanc, J., Mikulis, D. J., and Fisher, J. A. (2015). The dynamics of cerebrovascular reactivity shown with transfer function analysis. Neuroimage 114, 207–216. doi: 10. 1016/j. neuroimage. 2015. 04. 029

Duvernoy, H. M., Delon, S., and Vannson, J. L. (1981). Cortical blood vessels of the human brain. Brain Res. Bull. 7, 519–579. doi: 10. 1016/0361-9230(81)90007-1

Edvinsson, L., MacKenzie, E. T., and McCulloch, J. (2002). Cerebral Blood Flow and Metabolism. New York, NY: Raven Press.

Ellis, M. J., Ryner, L. N., Sobczyk, O., Fierstra, J., Mikulis, D. J., Fisher, J. A., et al. (2016). Neuroimaging assessment of cerebrovascular reactivity in concussion: current concepts, methodological considerations, and review of the literature. Front. Neurol. 7: 61. doi: 10. 3389/fneur. 2016. 00061

Fierstra, J., Poublanc, J., Han, J. S., Silver, F., Tymianski, M., Crawley, A. P., et al. (2010). Steal physiology is spatially associated with cortical thinning. J. Neurol. Neurosurg. Psychiatry 81, 290–293. doi: 10. 1136/jnnp. 2009. 188078

Fisher, J. A. (2016). The CO2 stimulus for cerebrovascular reactivity: fixing inspired concentrations vs. targeting end-tidal partial pressures. J. Cereb. Blood Flow Metab. 36, 1004–1011. doi: 10. 1177/0271678X16639326

Fisher, J. A., Venkatraghavan, L., and Mikulis, D. J. (2018). Magnetic resonance imaging–based cerebrovascular reactivity and hemodynamic reserve. Stroke 49, 2011–2018. doi: 10. 1161/STROKEAHA. 118. 021012

Frederick, B., Nickerson, L. D., and Tong, Y. (2012). Physiological denoising of BOLD fMRI data using regressor interpolation at progressive time delays (RIPTiDe) processing of concurrent fMRI and near-infrared spectroscopy (NIRS). Neuroimage 60, 1913–1923. doi: 10. 1016/j. neuroimage. 2012. 01. 140

Geranmayeh, F., Wise, R. J. S., Leech, R., and Murphy, K. (2015). Measuring vascular reactivity with breath-holds after stroke: a method to aid interpretation of group-level BOLD signal changes in longitudinal fMRI studies. Hum. Brain Mapp. 36, 1755–1771. doi: 10. 1002/hbm. 22735

Iadecola, C., and Nedergaard, M. (2007). Glial regulation of the cerebral microvasculature. Nat. Neurosci. 10, 1369–1376. doi: 10. 1038/nn2003

Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W., and Smith, S. M. (2012). Fsl. Neuroimage 62, 782–790. doi: 10. 1016/j. neuroimage. 2011. 09. 015

Juttukonda, M. R., and Donahue, M. J. (2019). Neuroimaging of vascular reserve in patients with cerebrovascular diseases. Neuroimage 187, 192–208. doi: 10. 1016/j. neuroimage. 2017. 10. 015

Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. P. W. (2010). Elastix: a toolbox for intensity-based medical image registration. IEEE Trans. Med. Imaging 29, 196–205. doi: 10. 1109/TMI. 2009. 2035616

Len, T. K., Neary, J. P., Asmundson, G. J. G., Goodman, D. G., Bjornson, B., and Bhambhani, Y. N. (2011). Cerebrovascular reactivity impairment after sport-induced concussion. Med. Sci. Sports Exerc. 43, 2241–2248. doi: 10. 1249/MSS. 0b013e3182249539

Liu, P., Welch, B. G., Li, Y., Gu, H., King, D., Yang, Y., et al. (2017). Multiparametric imaging of brain hemodynamics and function using gas-inhalation MRI. Neuroimage 146, 715–723. doi: 10. 1016/j. neuroimage. 2016. 09. 063

Liu, P. B., De Vis, J., and Lu, H. (2018). Cerebrovascular reactivity (CVR) MRI with CO2 challenge: a technical review. Neuroimage 187, 104–115. doi: 10. 1016/j. neuroimage. 2018. 03. 047

Lu, H., Golay, X., Pekar, J. J., and Van Zijl, P. C. M. (2003). Functional magnetic resonance imaging based on changes in vascular space occupancy. Magn. Reson. Med. 50, 263–274. doi: 10. 1002/mrm. 10519

Mandell, D. M., Han, J. S., Poublanc, J., Crawley, A. P., Stainsby, J. A., Fisher, J. A., et al. (2008a). Mapping cerebrovascular reactivity using blood oxygen level-dependent MRI in patients with arterial steno-occlusive disease: comparison with arterial spin labeling MRI. Stroke 39, 2021–2028. doi: 10. 1161/STROKEAHA. 107. 506709

Mandell, D. M., Han, J. S., Poublanc, J., Crawley, A. P., Kassner, A., Fisher, J. A., et al. (2008b). Selective reduction of blood flow to white matter during hypercapnia corresponds with leukoaraiosis. Stroke 39, 1993–1998. doi: 10. 1161/STROKEAHA. 107. 501692

Mark, C. I., and Pike, G. B. (2012). Indication of BOLD-specific venous flow-volume changes from precisely controlled hyperoxic vs. hypercapnic calibration . J. Cereb. Blood Flow Metab. 32, 709–719. doi: 10. 1038/jcbfm. 2011. 174

Mark, C. I., Slessarev, M., Ito, S., Han, J., Fisher, J. A., and Pike, G. B. (2010). Precise control of end-tidal carbon dioxide and oxygen improves BOLD and ASL cerebrovascular reactivity measures. Magn. Reson. Med. 64, 749–756. doi: 10. 1002/mrm. 22405

McDonald, D. (1974). Blood Flow in Arteries , 2nd Edn. London: Edward Arnold.

Moody, D. M., Bell, M. A., and Challa, V. R. (1990). Features of the cerebral vascular pattern that predict vulnerability to perfusion or oxygenation deficiency: an anatomic study. Am. J. Neuroradiol. 11, 431–439.

Ogawa, S., Lee, T. M., Kay, A. R., and Tank, D. W. (1990). Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc. Natl. Acad. Sci. U. S. A. 87, 9868–9872. doi: 10. 1073/pnas. 87. 24. 9868

Ogawa, S., Menon, R. S., Tank, D. W., Kim, S. G., Merkle, H., Ellermann, J. M., et al. (1993). Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model. Biophys. J. 64, 803–812. doi: 10. 1016/S0006-3495(93)81441-3

Paulson, O. B., Strandgaard, S., and Edvinsson, L. (1990). Cerebral autoregulation. Cerebrovasc. Brain Metab. Rev. 2, 161–192.

Poublanc, J., Crawley, A. P., Sobczyk, O., Montandon, G., Sam, K., Mandell, D. M., et al. (2015). Measuring cerebrovascular reactivity: the dynamic response to a step hypercapnic stimulus. J. Cereb. Blood Flow Metab. 35, 1746–1756. doi: 10. 1038/jcbfm. 2015. 114

Ringelstein, E. B., Sievers, C., Ecker, S., Schneider, P. A., and Otis, S. M. (1988). Noninvasive assessment of CO2-induced cerebral vasomotor response in normal individuals and patients with internal carotid artery occlusions. Stroke 19, 963–969. doi: 10. 1161/01. STR. 19. 8. 963

Rostrup, E., Law, I., Blinkenberg, M., Larsson, H. B., Born, A. P., Holm, S., et al. (2000). Regional differences in the CBF and BOLD responses to hypercapnia: a combined PET and fMRI study. Neuroimage 11, 87–97. doi: 10. 1006/nimg. 1999. 0526

Sam, K., Crawley, A. P., Conklin, J., Poublanc, J., Sobczyk, O., Mandell, D. M., et al. (2016). Development of white matter hyperintensity is preceded by reduced cerebrovascular reactivity. Ann. Neurol. 80, 277–285. doi: 10. 1002/ana. 24712

Shamonin, D. P., Bron, E. E., Lelieveldt, B. P. F., Smits, M., Klein, S., and Staring, M. (2014). Fast parallel image registration on CPU and GPU for diagnostic classification of Alzheimer’s disease. Front. Neuroinform. 7: 50. doi: 10. 3389/fninf. 2013. 00050

Slessarev, M., Han, J., Mardimae, A., Prisman, E., Preiss, D., Volgyesi, G., et al. (2007). Prospective targeting and control of end-tidal CO2 and O2 concentrations. J. Physiol. 581, 1207–1219. doi: 10. 1113/jphysiol. 2007. 129395

Sobczyk, O., Battisti-Charbonney, A., Fierstra, J., Mandell, D. M., Poublanc, J., Crawley, A. P., et al. (2014). A conceptual model for CO2-induced redistribution of cerebral blood flow with experimental confirmation using BOLD MRI. Neuroimage 92, 56–68. doi: 10. 1016/j. neuroimage. 2014. 01. 051

Sobczyk, O., Battisti-Charbonney, A., Poublanc, J., Crawley, A. P., Sam, K., Fierstra, J., et al. (2015). Assessing cerebrovascular reactivity abnormality by comparison to a reference atlas. J. Cereb. Blood Flow Metab. 35, 213–220. doi: 10. 1038/jcbfm. 2014. 184

Thomas, B. P., Liu, P., Park, D. C., van Osch, M. J. P., and Lu, H. (2014). Cerebrovascular reactivity in the brain white matter: magnitude, temporal characteristics, and age effects. J. Cereb. Blood Flow Metab. 34, 242–247. doi: 10. 1038/jcbfm. 2013. 194

Tong, Y., and Frederick, B. B. (2012). Concurrent fNIRS and fMRI processing allows independent visualization of the propagation of pressure waves and bulk blood flow in the cerebral vasculature. Neuroimage 61, 1419–1427. doi: 10. 1016/j. neuroimage. 2012. 03. 009

Tripp, B. C., Smith, K., and Ferry, J. G. (2001). Carbonic anhydrase: new insights for an ancient enzyme. J. Biol. Chem. 276, 48615–48618. doi: 10. 1074/jbc. R100045200

van der Zande, F. H., Hofman, P. A., and Backes, W. H. (2005). Mapping hypercapnia-induced cerebrovascular reactivity using BOLD MRI. Neuroradiology 47, 114–120. doi: 10. 1007/s00234-004-1274-3

van Niftrik, C. H. B., Piccirelli, M., Bozinov, O., Pangalu, A., Fisher, J. A., Valavanis, A., et al. (2017). Iterative analysis of cerebrovascular reactivity dynamic response by temporal decomposition. Brain Behav. 7: e00705. doi: 10. 1002/brb3. 705

Xu, F., Liu, P., Pascual, J. M., Xiao, G., and Lu, H. (2012). Effect of hypoxia and hyperoxia on cerebral blood flow, blood oxygenation, and oxidative metabolism. J. Cereb. Blood Flow Metab. 32, 1909–1918. doi: 10. 1038/jcbfm. 2012. 93

Yoon, S. H., Zuccarello, M., and Rapoport, R. M. (2000). Reversal of hypercapnia induces endothelin-dependent constriction of basilar artery in rabbits with acute metabolic alkalosis. Gen. Pharmacol. Vasc. Syst. 35, 333–340. doi: 10. 1016/S0306-3623(02)00112-X

Zhang, Y., Brady, M., and Smith, S. (2001). Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans Med Imaging 20, 45–57. doi: 10. 1109/42. 906424