Rapid cerebrovascular reactivity mapping: Enabling vascular reactivity information to be routinely acquired

Abstract Cerebrovascular reactivity mapping (CVR), using magnetic resonance imaging (MRI) and carbon dioxide as a stimulus, provides useful information on how cerebral blood vessels react under stress. This information has proven to be useful in the study of vascular disorders, dementia and healthy ageing. However, clinical adoption of this form of CVR mapping has been hindered by relatively long scan durations of 7–12 min. By replacing the conventional block presentation of carbon dioxide enriched air with a sinusoidally modulated stimulus, the aim of this study was to investigate whether more clinically acceptable scan durations are possible. Firstly, the conventional block protocol was compared with a sinusoidal protocol of the same duration of 7 min. Estimates of the magnitude of the CVR signal (CVR magnitude) were found to be in good agreement between the stimulus protocols, but estimates of the relative timing of the CVR response (CVR phase) were not. Secondly, data from the sinusoidal protocol was reanalysed using decreasing amounts of data in the range 1–6 min. The CVR magnitude was found to tolerate this reduction in scan duration better than CVR phase. However, these analyses indicate that scan durations in the range of 3–5 min produce robust data. HighlightsCerebrovascular reactivity mapping is a promising way to stress test the brain.Typically imaging times are relatively long compared with complementary methods.In this study we demonstrate a rapid cerebrovascular reactivity mapping technique.This is achieved by using a sinusoidally modulated hypercapnia challenge.This is shown to give equivalent results to the conventional block paradigm.


Introduction 43 44
Whilst there are many clinical techniques to investigate cerebral physiology at rest, 45 methods to examine the brain during dynamic changes in demand are few in number. 46 Resting perfusion measurement techniques using Computed Tomography (CT) or 47 Magnetic Resonance Imaging (MRI) are able to confirm that the tissue is well supplied 48 at rest, but they do not provide information on how the vasculature will react during 49 increases in demand. By introducing a physiological challenge this aspect of 50 cerebrovascular health can be revealed. There are two main vasoactive challenges that 51 can be used to provide a stress test of the cerebral vasculature: acetazolamide and 52 carbon dioxide. Acetazolamide provides a large physiological challenge by maximally 53 dilating cerebral arterioles. However, its effects are relatively long lasting (~30 mins) 54 and it therefore has the potential to interfere with the accuracy of other diagnostic 55 measurements made in the same session. In contrast, carbon dioxide may have a weaker 56 vasodilatory effect, but it is rapidly cleared by the lungs enabling multiple short 57 repeated challenges to be performed without interfering with subsequent measurements. 58 When combined with an imaging modality this carbon dioxide stimulus, often termed a 59 hypercapnia challenge, enables cerebrovascular reactivity (CVR) to be mapped at the 60 tissue level. The most common imaging modality for this purpose is MRI. 61 Measurements have been performed using both Arterial Spin Labelling (ASL) 1 and 62 Blood Oxygenation Level Dependent (BOLD) 2 weighted techniques, each with their 63 own advantages and disadvantages. BOLD weighted imaging offers higher spatial and 64 temporal resolution and a higher signal to noise ratio (SNR) than ASL, whereas ASL 65 provides quantitative measurements of perfusion change and BOLD does not. In 66 practice BOLD weighted measurements are more commonly used due to their higher 67 sensitivity, more widespread availability and ease of analysis. 68 The information provided by CVR measurements has proven useful in a range 69 of applications. CVR mapping has perhaps made the biggest impact in the study of 70 vascular disorders of the brain. In patients with unilateral carotid stenosis it has been 71 used to show that CVR is reduced in both the affected and unaffected hemisphere 3 . In 72 Moyamoya patients, who share similar pathological narrowing of cerebral vessels, 73 BOLD based CVR has been shown to give comparable results with ASL based CVR 1 . 74 CVR measurements have also been acquired in difficult patient populations such as in 75 subarachnoid haemorrhage 4 or in rare conditions such as MELAS 5 where in both cases 76 only limited biomarkers are available. CVR has also shown promise as a biomarker in 77 dementia and healthy ageing research. In Alzheimer's disease CVR has been shown to 78 be altered in patients with different variants of the APOE gene 6 . CVR has also been 79 used to investigate the effect of healthy ageing on brain haemodynamics 7 . 80 However, the clinical adoption of hypercapnia based CVR mapping has so far 81 been hindered by two main challenges: (i) relatively long scan times and (ii) long 82 patient preparation times. Clinical MRI examinations are typically split into discrete 83 diagnostic measurements of approximately 5 minutes in duration. However, currently 84 CVR examinations are between 7 and 12 minutes in duration, meaning that 85 compromises in diagnostic protocols must be made to accommodate such 86 measurements. Preparation of patients adds an additional time penalty that is dependent 87 on the means of gas administration and/or gas sampling. However, in this study we 88 concentrate on the question of reducing the duration of the CVR examination. 89 Typically, hypercapnic gas challenges are presented in blocks interleaved with a 90 baseline air breathing condition. One popular gas protocol used by researchers in 91 Toronto involves the use of two hypercapnic blocks of differing duration; an initial 92 short block followed by a second much longer block 8 . This protocol has also been used 93 to derive information about the relative timing of the CVR response using Transfer 94 Function Analysis (TFA) 9 . This timing information has been suggested to be reflective 95 of regional variations in blood arrival time 10 or in the rate at which the CVR response 96 develops 9 . However, this block paradigm does not lend itself to easy scan time 97 reduction as it typically relies on long hypercapnic blocks to minimise the effect of 98 regional differences in vascular delays. Alternatively, in the limit of short blocks a 99 sinusoidally varying stimulus has several useful properties 10  prospective algorithm for the targeting of blood gases 12 . Subjects were fitted with a 148 sequential gas delivery (SGD) breathing circuit, which was sealed to the face using 149 adhesive tape (Tegaderm, 3M Healthcare, St. Paul, Minnesota, USA) to prevent the 150 entrainment of room air. Subjects were then asked to sit comfortably outside the scanner 151 whilst the prospective algorithm was calibrated. The subject's specific resting end-tidal 152 carbon dioxide and oxygen partial pressures (PETCO 2 and PETO 2 , respectively) were 153 first determined by averaging the first ten or so breaths into the breathing mask. The 154 calibration consisted of initiating a baseline targeted at the subject's specific resting 155 PETCO 2 and PETO 2 . Initial estimates of the subject's VCO 2 and VO 2 (their resting 156 production/consumption of CO 2 /O 2 ) based on height, weight and age were refined until 157 a constant baseline was achieved with minimal drift. A brief presentation of the 158 sinusoidal hypercapnic stimulus was then administered to the subject to prepare them 159 for the sensations they might encounter whilst in the scanner. Gas stimulus protocols (Fig. 1a) were tailored to each subject's resting PETCO 2 161 and PETO 2 baselines. Modulations in PETCO 2 were targeted relative to baseline, whilst 162 PETO 2 was targeted to remain constant at baseline throughout. The Toronto gas protocol 163 consisted of two blocks of hypercapnia; one of 45 s duration and the other 120 s. The During the presentation of each 7 minute stimulus protocol, imaging was 173 performed using the 7 minute BOLD-weighted EPI pulse sequence described above. 174 Synchronisation of the respiratory and imaging protocols was manually initiated. 175 176

Image Analysis 177
The basic image analysis pipeline is summarised in Fig. 1b. Imaging data from the two 178 stimulus protocols were first submitted to a common pre-processing stage including 179 motion correction 14 , fieldmap based EPI unwarping 15 , slice timing correction, removal 180 of non-brain tissue 16  changes in the brain, due to intersubject differences in blood arrival time. This was 194 achieved by extracting the mean whole brain BOLD signal change, using the mask 195 generated during the brain extraction process, and finding the maximum cross-196 correlation with the PETCO 2 data interpolated to the BOLD sampling frequency. High 197 pass temporal filtering with the same cut-off as the imaging data was then applied and 198 the temporal derivative included as an additional regressor. To provide CVR phase estimates that are comparable with the Toronto protocol, maps 212 were produced using the whole brain phase as a reference i.e. equivalent to the 213 synchronisation performed for the Toronto protocol. This was achieved by extracting a 214 whole brain BOLD signal timecourse and estimating the CVR phase using Eq. (3). This 215 global phase value was then subtracted from the CVR phase estimates from each voxel. 216 For the Toronto protocol Transfer Function Analysis (TFA) 9 was used to 217 estimate CVR phase. Since this is a Fourier based method, and as such no standard tools The standard deviation in the CVR magnitude ( :NO P QRS ) can be calculated using Eq. 230 (7) However, the absolute error in CVR magnitude does not take into account the large 237 differences in the scale of CVR magnitude particularly between grey and white matter. 238 Under these conditions it is preferable to represent the standard deviation as a fraction 239 of CVR magnitude by calculating the relative standard deviation (RSD).  The normocapnic PETCO 2 (PETCO 2,norm ) and the peak change in PETCO 2 (DPETCO 2 ) 248 were measured for each protocol and subject. First, the automated end-tidal picking 249 routine was manually inspected and corrected where necessary, using software provided 250 with the RespirAct, to ensure accurate PETCO 2 estimation. Linear regression was then 251 performed using a prototypical model of the stimulus (i.e. a boxcar temporally aligned 252 with the PETCO 2 data or sine/cosine pair at 16.7 mHz) and a constant term; the former 253 representing DPETCO 2 and the latter PETCO 2,norm . 254 255

Comparison between Toronto and Sinusoid protocols 256
Investigations were made to compare the sensitivity and equivalence of the two 257 protocols. Firstly, the number of statistically significant voxels that passed the corrected 258 voxel level threshold (p<0.05) were counted for each subject and protocol. In order to 259 control for potential differences in subject motion between protocols the mean motion 260 as calculated during the motion correction step was extracted 14 . Mean motion represents 261 the mean-root-mean square displacement between adjacent time points and has been 262 shown to be associated with reduced temporal signal to noise ratio (SNR) in fMRI 263 experiments 22 . It is important that the effective temporal SNR is fairly constant across 264 protocols in order to isolate the effect of different stimulus paradigm rather than 265 difference in noise characteristics. 266 Next the equivalence of the CVR estimates produced by the Toronto and 267 Sinusoid protocols was explored. Cortical ROIs were taken from the co-registered MNI 268 atlas and further refined using a grey matter mask (grey matter segmentation partial 269 volume estimate thresholded at 0.5). The mean and standard deviation of CVR 270 magnitude and phase values extracted from these ROIs were calculated for each subject. The effect of reducing the amount of data used to measure the CVR parameters using 276 the Sinusoid protocol was explored. This was achieved by truncating the data set at scan 277 durations between 1 minute (30 volumes) and 6 minutes (150 volumes) in steps of 1 278 minute (30 volumes). The full analysis, including pre-processing and GLM analysis, 279 was performed for each truncated dataset, and maps of CVR magnitude and phase were 280 generated. High pass temporal filtering at 100 s was disabled for the 1 minute scan 281

duration. 282
The impact of scan time reduction was investigated at the level of whole brain 283 grey matter. An ROI was generated based on voxels in the Toronto protocol that met the 284 statistical threshold and was then refined using a grey matter mask (grey matter 285 segmentation partial volume estimate thresholded at 0.5). The mean and standard 286 deviation of the CVR magnitude and phase values in this ROI were calculated for each 287 scan duration and subject. In addition, the number of statistically significant voxels 288 from the GLM analysis were calculated for each scan duration and subject. Although within the accuracy of the end-tidal gas targeting system, on average baseline 302 PETCO 2 was 0.9 mmHg higher and DPETCO 2 was 0.7 mmHg lower for the Toronto 303 protocol compared with the Sinusoid protocol. Paired two tailed T-tests showed that 304 these differences were significant at p<0.001 and p<0.01, respectively. Although the 305 target of 10 mmHg PETCO 2 change was not attained, all subjects experienced similar 306 PETCO 2 changes within a standard deviation of 1.2 mmHg. 307 In order to compare the sensitivity of the two protocols the number of voxels 308 that passed the voxel level statistical threshold were counted. Table 2   Sinusoid protocol compared with the Toronto protocol (Fig. 2a/e). However, similar 325 grey matter to white matter contrast is observed across CVR magnitude (Fig. 2b/f)  with the amount of data included in the analysis (Fig. 4a-c). CVR magnitude maps look 339 qualitatively very similar across scan durations (Fig. 4d-f), whilst CVR phase maps 340 show some small differences when only 3 minutes of data are used (Fig. 4g-i). This is 341 further explored across all subjects in Fig. 5. The group mean of the mean CVR 342 magnitude across the ROI (Fig. 5a) was seen to increase slightly with reduced scan 343 time. A two-way ANOVA test showed that not all group means were equal (p<0.001), 344 but a paired multiple comparisons analysis showed that only the 1 minute duration was 345 significantly different (p<0.001) to the full data set. A similar increase was observed for 346 the group mean standard deviation of the CVR magnitude (Fig. 5b). Again group means 347 were not all equal (p<0.001) and only the 1 minute duration was significantly different 348 (p<0.001) to the full data set. The group mean of the mean CVR phase displayed little 349 variation (Fig. 5c) and was not significant (p=0.99). The group mean of the standard 350 deviation of the CVR phase showed the largest increase with time reduction (Fig. 5d). CVR mapping using MRI is a promising clinical tool for stress testing cerebral blood 363 vessels. However, there are several aspects of current CVR mapping techniques that 364 limit their clinical application. In this study we concentrated on minimising the scan 365 duration in order to improve the compatibility of this technique with clinical workflows. 366 By utilising a sinusoidally modulated hypercapnia challenge we were able to 367 retrospectively examine the effect of reducing the scan duration. In turn we were able to 368 demonstrate that good quality CVR maps can be acquired in around 3 minutes, rather 369 than the more typical 7 minutes. Furthermore, we demonstrated that CVR maps 370 generated using this sinusoidal stimulus are equivalent to the more conventional block 371 paradigm. 372 373

Stimulus performance 374
Gas stimulus protocols were tailored to each subject's individual baseline PETCO 2 . This 375 baseline was maintained within 1 mmHg between the protocols, but was statistically 376 significantly different none the less. Since the vascular response to carbon dioxide is 377 sigmoidal in form 23 it is possible that the two protocols were performed on different 378 parts of the sigmoid curve. However, in practice maintaining baseline PETCO 2 to a 379 higher accuracy than 1 mmHg is challenging. Both stimulus protocols elicited PETCO 2 380 changes of the order of 7 mmHg, which should place the vascular response within the 381 linear regime of the sigmoidal vascular response 23 . A statistically significant difference 382 in DPETCO 2 was observed, although well within the accuracy of the gas delivery system 383 at 0.7 mmHg. Unfortunately, subject tolerance of both stimulus protocols was not 384 recorded. However, none of the subjects reported discomfort from either stimulus 385 protocol. 386 387

Comparison between Toronto and Sinusoid protocols 388
Before discussing the results of the protocol comparison it is useful to consider the 389 advantages and disadvantages of the two stimulus protocols. As a more conventional 390 block paradigm, the Toronto protocol is an easier stimulus to implement. However, 391 synchronisation of the measured PETCO 2 experienced by the subject with the BOLD 392 data is critical to the accurate estimation of magnitude of the CVR change. Whilst this is 393 practical when synchronising with the whole brain average BOLD signal, the presence 394 of local variations in vascular delays opens up the possibility of the underestimation of 395 CVR magnitude, or the measurement of spurious negative signal changes 9 . This can be 396 mitigated by the use of long blocks of hypercapnia, but therefore limits how short the 397 total acquisition can be. In contrast sinusoidal stimuli have seen only limited 398 application 24 . To be clinically practical, computer controlled gas delivery is desirable, 399 although a manual implementation has also been demonstrated 25 . In the previous 400 implementation a Fourier analysis approach was used based on in-house analysis 401 software 10 . However, in this study we reframe the analysis in order to use free and 402 widely available fMRI analysis tools 15 . One of the inherent advantages of the Sinusoid 403 protocol is its insensitivity to local vascular delays. This is due to the fact that any 404 sinusoid with arbitrary amplitude and phase can be represented by the sum of a sine and 405 a cosine with different amplitudes. In this context only the frequency of the stimulus is 406 important, rather than the timing. This enables the short time-period stimuli used in this 407 work, which would otherwise result in underestimates, or artefactual negative values, of 408 CVR magnitude. Finally, simple sine and cosine functions were used as regressors in 409 the GLM analysis, rather than the measured PETCO 2 timecourse. Despite this, the 410 number of statistically significant voxels (Table 2) was not significantly different to the 411 Toronto protocol across the group. This approach simplifies the analysis as the PETCO 2 412 data does not need to be transferred to the analysis computer and synchronised with the 413 BOLD data. 414 Maps of CVR magnitude were found to have qualitatively similar features (Fig.  415 2b/f). Furthermore, quantitative analysis comparing regional estimates of CVR 416 magnitude at the level of the cerebral lobes demonstrated significant correlation 417 between protocols (Fig. 3a) and equivalence of the absolute values (Fig. 3c). By 418 utilising estimates of the variance in the parameter estimates and the residuals from the 419 GLM analysis it was possible to evaluate the standard deviation of the CVR magnitude 420 measurement. By presenting this as the relative standard deviation (RSD), with respect 421 to the mean, the large degree of heterogeneity in CVR magnitude is controlled, 422 revealing more noisy estimates in white matter (Fig. 2c/g). However, the standard 423 deviation in CVR magnitude, as estimated using Eq. (6), should be viewed with caution 424 as it relies on the variance in S mean from the residuals. The residuals include 425 contributions from motion and physiological noise that aren't included in the GLM 426 model. Therefore, large amounts of motion will result in larger values of RSD, in effect 427 being categorised as noise. 428 Maps of CVR phase also showed clear similarities between methods (Fig. 2d/h). 429 Quantitative analysis showed that the protocols were significantly correlated (Fig. 3b). 430 However, the slope between the protocols varied in the range -0.63 to 1.44 and was 431 trending towards being significantly different to unity (the value expected for 432 equivalence). This disparity may be due to the properties of phase, particularly the 433 incorrect unwrapping of phase values close to ±2p. It is also possible that CVR phase 434 behaves differently between the protocols given the large difference in stimulus 435 duration between the Toronto (45-120 s above baseline PETCO 2 ) and the Sinusoid 436 protocol (30 s above baseline PETCO 2 in each 60 s cycle). Maps of the standard 437 deviation of the CVR phase measurement were estimated for the Sinusoid protocol by 438 propagating the errors in the parameter estimates (Fig. 2i). RSD was not used because 439 of the relative nature of phase i.e. the reference phase is arbitrary. It is not possible to 440 estimate the standard deviation of the Toronto CVR phase due to the use of TFA, 441 although an analysis of the effects of noise was included in the original paper 9 . duration, there appears to be relatively little effect on maps of CVR magnitude at scan 448 durations of 3 or 5 minutes, and only small differences in the 3 minute CVR phase map 449 (Fig. 4). Quantitative analysis confirmed that the number of above threshold voxels in 450 the Z-statistic maps was significantly different in the 1, 2 and 3 minute scan durations 451 when compared with the full data set. Furthermore, the results in Fig. 5e are consistent  452 with the expected reduction in sensitivity demonstrated for fMRI analysis 26 . CVR 453 magnitude was only minimally affected, consistent with the qualitative assessment 454 above, with only the 1 minute scan duration found to be significantly different to the 455 full data set for both the mean and standard deviation across the ROI. Estimates of CVR 456 phase appear to be more strongly affected by scan time reduction. Whilst the mean over 457 the ROI was not significantly different as a function of scan duration, the standard 458 deviation was significantly different for the 1 and 2 minute durations. This increased 459 sensitivity of CVR phase to scan time reduction, relative to CVR magnitude, was 460 highlighted by normalising Fig. 5b and Fig. 5d to their 7 minute values. The standard 461 deviation of CVR phase over the ROI clearly increases much more rapidly for CVR 462 phase than CVR magnitude, necessitating longer scan duration where CVR phase is a 463 priority. In summary, the results would suggest that scan time could be reduced to be in 464 the range of 3 to 5 minutes with minimal effect on the resulting maps of CVR. baseline with minimal drift using the prospective end-tidal gas targeting system used in 471 this study, the system must be calibrated to each individual. In our hands this process 472 takes about 30 minutes, including application of the mask, although it may be shorter 473 for more experienced users. However, there is a great deal of scope for time reduction 474 by automating this calibration step, which will likely be implemented in future hardware 475 developments. In the short term, using the current hardware, it should still be possible to 476 realise shorter patient preparation times, albeit at the expense of a small amount of 477 PETCO 2 baseline drift. As noted above, the gas targeting system provides an initial setup 478 based on the height, weight and age of the subject. The calibration process enables this 479 setup to be refined. However, the initial estimate results in only small amounts of drift 480 in many cases. Small amounts of drift are removed by the high-pass temporal filter and 481 will not affect estimates of CVR so long as the vasculature remains in the linear regime