Approaches to Brain Stress Testing: BOLD Magnetic Resonance Imaging with Computer-Controlled Delivery of Carbon Dioxide

Background An impaired vascular response in the brain regionally may indicate reduced vascular reserve and vulnerability to ischemic injury. Changing the carbon dioxide (CO2) tension in arterial blood is commonly used as a cerebral vasoactive stimulus to assess the cerebral vascular response, changing cerebral blood flow (CBF) by up to 5–11 percent/mmHg in normal adults. Here we describe two approaches to generating the CO2 challenge using a computer-controlled gas blender to administer: i) a square wave change in CO2 and, ii) a ramp stimulus, consisting of a continuously graded change in CO2 over a range. Responses were assessed regionally by blood oxygen level dependent (BOLD) magnetic resonance imaging (MRI). Methodology/Principal Findings We studied 8 patients with known cerebrovascular disease (carotid stenosis or occlusion) and 2 healthy subjects. The square wave stimulus was used to study the dynamics of the vascular response, while the ramp stimulus assessed the steady-state response to CO2. Cerebrovascular reactivity (CVR) maps were registered by color coding and overlaid on the anatomical scans generated with 3 Tesla MRI to assess the corresponding BOLD signal change/mmHg change in CO2, voxel-by-voxel. Using a fractal temporal approach, detrended fluctuation analysis (DFA) maps of the processed raw BOLD signal per voxel over the same CO2 range were generated. Regions of BOLD signal decrease with increased CO2 (coded blue) were seen in all of these high-risk patients, indicating regions of impaired CVR. All patients also demonstrated regions of altered signal structure on DFA maps (Hurst exponents less than 0.5; coded blue) indicative of anti-persistent noise. While ‘blue’ CVR maps remained essentially stable over the time of analysis, ‘blue’ DFA maps improved. Conclusions/Significance This combined dual stimulus and dual analysis approach may be complementary in identifying vulnerable brain regions and thus constitute a regional as well as global brain stress test.


Introduction
The brain is exquisitely sensitive to alterations in carbon dioxide (CO 2 ) tension in the blood. [1] Much time and effort is spent in neuroanesthesia and neurocritical care to control arterial CO 2 tension because a higher CO 2 increases intracranial pressure and volume. Increased CO 2 in normal circumstances increases cerebral blood flow (CBF) by up to 5-11percent/mmHg [2] with a concomitant increase in cerebral blood volume. Alterations in this tight relationship between CO 2 tension and CBF occur with brain injury and disease states. Frequently cerebrovascular reactivity (CVR) to CO 2 is diminished or in certain circumstances actually reversed -a phenomenon termed 'cerebral steal'. [3] [4] Patients with these alterations are at far greater risk of a poor outcome following surgical intervention; for example, stroke. Thus, manipulation of CO 2 to prospectively assess patient risk of cerebral injury has been used for some time. These include increases in CO 2 by breath-holding or administration of acetazolamide to increase intracellular CO 2 and hydrogen ion concentration. [5] To date, blood velocity changes following CO 2 challenge are routinely assessed by transcranial Doppler. [6] Improving the reproducibility and precision of the CO 2 challenge for cerebral stress testing married to magnetic resonance imaging (MRI) for true regional assessment represents a significant step forward. [7] This objective, is in part, addressed in this study.
Higher field strength MRI using blood oxygen level dependent (BOLD) contrast has provided improved resolution in following changes in CBF and tissue oxygenation. Resting state [8] [9] and default mode [10] studies have helped define consciousness, and task-related BOLD imaging has demonstrated brain regional interconnectivity. [11] [12] Here we describe dynamic brain imaging in health and disease using carefully controlled alterations in CO 2 to help understand CVR. We have used two patterns of changes in CO 2 (a stepwise change in CO 2 and a ramp change) and analyzed the resultant MR-BOLD output as a voxel-based ratio of change in mean BOLD signal to change in CO 2 to determine regional CVR [13] -and a voxel-based fractal time analysis of the MR-BOLD signal using detrended fluctuation analysis (DFA). [14] The latter technique has been successfully applied to identify activated cortical areas with functional MRI. The similarities and differences between the two carefully controlled CO 2 challenges and analysis paradigms are outlined. Further understanding of these relationships should advance the concept of a 'brain stress test' to aid in diagnosis and treatment of brain disease.

Methods
All patients and volunteer subjects signed written informed consent with the various protocols approved by the Ethics Review Board of the University Health Network of the University of Toronto.

Patient Group
Eight patients with severe carotid stenosis were screened by questionnaire for demographic data and any contraindication to undergoing MRI. Patients had been referred to the CVR study Figure 1. Response to square wave CO 2 stimulus -CVR and DFA maps. These data are from a patient who underwent the square wave sequence and ramp sequence ( Figure 2) in the same sitting. Not shown in this image is the end-tidal O 2 tension which is stable over time at normal values (approximately 100 mmHg). The CO 2 stimulus is shown in the first column -start, middle and end during the two-minute stimulus. The first 45-second square wave pulse was to aid in time sequencing of the BOLD signal and CO 2 stimulus. Time is on the x-axis and end-tidal CO 2 in mmHg on the y-axis. The yellow and magenta vertical markers highlight the center of the CO 2 durations analyzed in the square wave sequence. The highlighted red lines are the points where the CO 2 tension was correlated to the MR-BOLD signal for the CVR analysis. The second column shows the corresponding CVR maps. The color key to the right of the image is upper red +0.56 arbitrary BOLD units/mmHg in CO 2 tension -deep blue 20.56. The third column shows the corresponding DFA maps. For the DFA maps the BOLD signal analysis was based on output as interpolated between the yellow and magenta markers. The color key to the right of the image is brown -1.5 Hurst units, pink -1.0 Hurst units, white -0.5 Hurst units and deep blue -0 Hurst units. 'Blue' DFA was defined as less than 0.5 Hurst units and indicates anti-persistent noise. The percentage of 'blue' DFA voxels noticeably decreased over time. The step change in end-tidal CO 2 is 8 mmHg. doi:10.1371/journal.pone.0047443.g001 group because of significant cerebrovascular pathology. After familiarization with the protocol and fit to a breathing circuit (see below) the patients were placed in the bore of a 3.0 Tesla MR unit (Signa; GE Healthcare, Milwaukee, WI) and imaged per protocol with a BOLD echo planar sequence to assess the regional brain response to alterations in CO 2 (see below).
Volunteer Subject Group -Two healthy volunteers were studied. After familiarization with the protocol they were fit to a breathing circuit (see below). They also had simultaneous bilateral recording of middle cerebral artery blood velocity by transcranial Doppler (Spencer Technologies, Seattle, WA), and bilateral frontal cerebral oximetry (Fore-Sight, CasMed, Branford, CT), continuous finger plethysmography (Nexfin, BMEYE, Amsterdam, The Netherlands), and continuous pulse oximetry (Nellcor, Covidien, Mansfield, MA). The hemodynamic responses to the changes in CO 2 were recorded. Response to ramp CO 2 stimulus -CVR and DFA maps. The CO 2 ramp stimulus is shown in the first column. The duration of the ramp stimulus increases by row. Not shown in this image is the end-tidal O 2 tension which is stable over time at normal values (approximately 100 mmHg). The yellow and magenta vertical markers highlight the center of the CO 2 durations analyzed in the square wave sequence. The highlighted red lines demonstrate the interpolated points where the CO 2 tension was correlated to the MR-BOLD signal for the DFA analysis. The CVR maps are shown in second column and the DFA maps in the third column. In the first time period CVR map (row 1, column 2) the deep saturation of red signal indicates initial very high responsiveness (steep slope) within these areas suggesting these voxels are on the steep linear portion of the sigmoidal CO 2 response curve. In the DFA map (row 1, column 3) these areas are noticeably brown with a Hurst exponent of 1.5 indicating very high time memory. The ramp change in end-tidal CO 2 examined is 8 mmHg. The color keys are as in Figure 1. doi:10.1371/journal.pone.0047443.g002 Model-based prospective end-tidal (MPET) gas breathing sequences The theory and application of the MPET approach (RespirAct -Thornhill Research Inc. Toronto, Canada) to deliver carefully controlled gas mixtures has been described in detail previously. [15] In brief, the basic function of the sequential rebreathing circuit (fit to the study subject) is to effectively limit the inspired gas to that of the output from the gas blender, independent of actual minute ventilation. This occurs because at any minute ventilation that exceeds the flow of gas from the blender the balance of inhaled gas is provided by previously exhaled gas from a circuit reservoir. The previously exhaled gas does not contribute to gas exchange as it has already equilibrated with the arterial blood. The flows and component gas concentrations of the gas entering the breathing circuit are calculated to target specific end-tidal CO 2 and O 2 tensions using algorithms described by Slessarev et al. [15] The timing and sequence of target end-tidal CO 2 and O 2 are prospectively entered into the MPET controller, which then implements them. With this system, the exhaled CO 2 values have been shown to precisely reflect those in the arterial blood, [16] the true independent variable for brain blood flow. This system is capable of generating highly reproducible square wave and ramp changes in CO 2 under constant O 2 tensions as used here (henceforth when discussing the gas mixtures programmed by the MPET we will note these as CO 2 or O 2 but imply end-tidal partial pressure of CO 2 or O 2 unless otherwise stated).

MRI sequences
MRI was performed on a 3.0-Tesla scanner with an 8-channel phased array head coil. T1-weighted anatomic images were acquired using a 3-dimensional spoiled gradient echo pulse sequence (whole brain coverage; matrix: 2566256; slice thickness: 2.2 mm; no interslice gap). BOLD MRI data were acquired with a T2*-weighted single-shot gradient echo pulse sequence with echo planar readout (field of view: 24624 cm; matrix: 64664; TR: 2000 ms; TE: 30 ms; flip angle: 85u; slice thickness: 5.0 mm; interslice gap: 2.0 mm; number of frames: 254). Post-hoc Analysis Patient Group: CVR maps Post-processing of the BOLD signal was per the usual techniques employing standard analysis of functional neuroimages (AFNI). [17] Custom-designed software time sequenced the CO 2 alterations to the BOLD output. Here CVR is defined as the change in BOLD signal from a mean baseline measurement to a mean measurement at the altered CO 2 /mean change in CO 2 tension between the two measurement periods. A full description is given in reference 13. The CVR was thresholded from 20.56 to +0.56 arbitrary BOLD units/mmHg change in CO 2 . Additionally, the raw MR output was imported into custom-written software (LabVIEW, National Instruments, Austin, TX). The CVR analysis was based on a two minute square wave alteration in CO 2 tension from 40-50 mmHg. 'Blue' CVR maps -a decrease in BOLD signal with the CO 2 stimulus (threshold 0 to 20.56) were quantified as the percent of total voxels for the axial slice. Three time periods were chosen at the start, middle and end of the sequence. The axial slice with the greatest number of voxels was chosen for analysis in each case. In the patient undergoing both square wave and ramp CO 2 protocols sequentially, the ramp data were analyzed in an incrementally increasing fashion over the same CO 2 range -see Figures 1 and 2 respectively.

Post-hoc Analysis Patient Group: DFA maps
DFA is a temporal fractal analysis approach. The raw BOLD signal sequenced to the CO 2 alterations noted above were analyzed using custom-written software (LabVIEW). The interpolated time periods of those used in the CVR maps were analyzed. DFA is a modified root mean square analysis of the magnitude of the noise signal fluctuations based on the granularity of the data (the time series is divided into boxes of equal length, with the longer box length customarily associated with greater signal fluctuation). A comprehensive description of the technique in relation to BOLD signal processing is given by Hu et al. [14] More generalized introductions to the technique are given by West [18] and Seely and Macklem. [19] A log-log plot of the relationship between the magnitude of fluctuation with increasing box length versus the box length is calculated. If the log-log transform of these data points have linear characteristics, the slope of this line is designated the Hurst exponent. The magnitude of the Hurst exponent is associated with well described noise patterns that have historically been identified with colors. Data were analyzed with a threshold for the Hurst exponent from 0 to 1.5. The color coding of the DFA maps were as follows -brown for a Hurst exponent in the range of 1.5 (Brownian noise), pink for a Hurst exponent of 1.0 (pink noise), white for a Hurst exponent of 0.5 (white noise) and blue for a Hurst exponent of 0 (anti-correlated noise). 'Blue' DFA maps -an anti-correlated time course of the BOLD signal with CO 2 stimulus (threshold 0.5 to 0.0) were quantified as percent of total voxels for the axial slice. Such 'blue' maps were compared and contrasted to the 'blue' CVR maps described above.

Post-hoc Analysis Volunteer Subject Group
The data were recorded and stored using a custom-designed digital acquisition system (LabVIEW). Data were analyzed following download to spreadsheet and time sequenced. The various hemodynamic responses to the same ramp sequence used in the MRI studies were analyzed. Heart rate, mean arterial pressure, middle cerebral artery velocity, and frontal cerebral oximetry responses were correlated to the change in CO 2 generated by the MPET gas blender.

Statistical Analysis
Data were analyzed using XLStat. Repeated measures ANOVA was used to compare CVR to DFA over the 3 time periods studied; p,0.05 considered significant by Tukey's test. Histogram analysis was fit to a maximal bin number of 25 for both CVR and DFA output. Correlations between hemodynamic variables versus CO 2 changes were assessed by linear regression analysis; p,0.05 considered significant.

Results
The demographics of the 8 patients with carotid stenosis are shown in Table 1. Angiographically documented severe cerebrovascular stenosis or occlusion was seen in all patients.
A comparison of CVR maps with DFA maps over time in one patient where both the square wave and ramp sequences were performed sequentially in one sitting are shown in Figures 1 and 2 Figures  for a description of differences). Visual inspection indicates that the 'blue' voxel counts resolve more rapidly over time with DFA. Also of note is the greater percentage of 'blue' voxels with the ramp protocol, especially at the outset, but, as in the square wave CO 2 challenge, resolution occurs in the later images. Also in the ramp sequence note the intense red saturation in the CVR map and the bleaching over time and the 'brown' saturation of high signal in the corresponding DFA map.
The density histograms collating the voxel-by-voxel output for each image as seen in Figures 1 and 2 Figure 5 shows results from a repeated measures ANOVA for data from the 8 patients imaged with steno-occlusive disease of their carotid arteries. A significant group 6 time interaction for percentage of 'blue' voxels over time is seen with p = 0.008. Within group comparisons reveal that the percentage of 'blue' voxels is significantly less with DFA compared to CVR at the end measurement period (p = 0.032 by Tukey's test). Figure 6 shows the hemodynamic effects on a 60 year old subject undergoing the same ramp sequence as used for the patient highlighted in Figure 2. A marked hyperdynamic blood pressure and cerebral blood flow velocity response is seen with increasing CO 2 that persists beyond the application of the CO 2 stimulus. By contrast, Figure 7 shows the hemodynamic effects in a 31year old subject with the ramp sequence; comparatively a much-dampened response is seen. Furthermore there is no increase in cerebral oxygen saturation to CO 2 stimulus.

Supporting Information -Movie S1
The accompanying movie shows the dynamic change in CVR calculated from the alteration in BOLD signal response to ramp changes in CO 2 tension over time.

Discussion
This study examines dynamic BOLD imaging of the brain and response to computer-controlled, repeatable CO 2 stimuli. The potential to quantify the cerebrovascular responsiveness and time course in patients at serious neurological risk from ischemic events can be appreciated. The use of two approaches to initiate a controlled change in CO 2 (square wave and ramp stimuli) and two analysis techniques of the BOLD signal (CVR and DFA) highlight how the brain regionally responds to changes in CO 2 that cannot be gleaned solely by one approach or one analysis technique. Our analysis of the dynamic response to the ramp stimulus indicates both similarities and differences to that seen with the more commonly reported square wave stimulus. We note that the square wave challenge provides a window into the dynamic effect of time on the cerebral vasculature response to a constant stimulus. Whereas, by contrast, the ramp challenge highlights the effect of increasing CO 2 alone as time is a controlled variable as the rate of change of the increasing CO 2 signal is constant. The dynamics of the response to CO 2 can be appreciated by examination of the first measurement period in the ramp sequence (row 1 in Figure 2). The intense red signal in the CVR map (middle column) and brown signal in the DFA map (right column) suggest that these voxels are demonstrating a CO 2 response centered about the steep linear portion of the sigmoidal cerebrovascular response curve at a Figure 6. Vigorous hemodynamic response to CO 2 ramp stimulus. Hemodynamic response to ramp CO 2 stimulus in a 60 year old healthy subject. Note the marked blood pressure and cerebral blood flow velocity response to increased CO 2 persisting past the stimulus peak. Legend abbreviations -etPCO 2 ; end-tidal partial pressure of carbon dioxide, MCAv; middle cerebral artery velocity, MAP; mean arterial pressure, HR; heart rate, right and left O 2 ; cerebral O 2 saturation. doi:10.1371/journal.pone.0047443.g006 CO 2 of 40 mmHg. On this portion of the CO 2 response curve cerebral blood flow velocity can change as much as 11percent/ mmHg increase in CO 2 tension. [2] The color fading seen in the latter two periods (rows 2 and 3) indicate a diminished responsiveness (a shallower slope) due to incorporation of BOLD signal averaging while on the plateau of the flow response to increases in CO 2 tension.
The generation of CVR maps using square wave CO 2 stimuli with these techniques have been previously reported. [20][21] Important therapeutic decisions have been based on longitudinal imaging in patients [7]; the reproducibility of the CO 2 stimulus allows such comparisons to be reliably made.
Detrended fluctuation analysis (DFA) is one approach to fractal signal processing. [18] [22] The technique has been applied successfully to functional MRI signal analysis. [14] The application of DFA to BOLD signal analysis is predicated on the premise that the technique is able to identify whether or not the signal examined has a 'time memory'. White noise has no memory -the signal is truly random -and has a Hurst exponent of 0.5 when the signal is processed by DFA. In contrast, many healthy physiological signals have a Hurst exponent greater than 0.5 and often close to 1.0 -pink noise. This is noise with a memory. The greater exponent indicates a repetitive pattern to the signal at every scale. [23] In the functional MRI study referenced activated voxels often demonstrated a Hurst exponent near unity. The interpretation is that the repetitive stimulus -say finger tapping -increases the BOLD signal in a reproducible manner during the course of the stimulus. The BOLD signal which is inherently noisy now has its noise altered by the superimposed stimulus -noise now with a memory of the activation initiated by the finger tapping. Examination of our data indicates that a reproducible CO 2 stimulus induces a robust noise memory -so much so that the Hurst exponent approached or exceeded 1.5 -Brownian noise. In a novel manner we also examined the distribution of 'blue' voxels with a Hurst exponent of less than 0.5. This exponent is defined as being associated with anti-correlated or anti-persistent noise and we interpret these findings as indicating a decreasing BOLD signal with increasing CO 2 -akin to the negative CVR slope indicating cerebral flow redistribution or steal depicted by 'blue' CVR voxels.
The DFA findings suggest that this technique also provides robust maps correlated to changes in BOLD signal intensity. The intensity of this signal is likely most influenced at the venular level. [24] The oxygen saturation of venous hemoglobin, venous volume, and metabolic rate can all influence the signal intensity, and recent work indicates that altered diffusion can also influence BOLD signal strength. [25] DFA demonstrates a dynamic response as opposed to CVR, which is a more static ratio of net signal change to the change in CO 2 . DFA is calculated from raw signal over time, while the signal for CVR is a meaned value over time at different CO 2 levels. As such, the changes in BOLD signal as measured by DFA may be tracking dynamic changes in oxygen extraction fraction from the blood. Consistent with this hypothesis is the finding of an increase in signal memory (increasing Hurst exponent) with activated voxels in functional MRI studies as a surrogate for increased venous oxygenation (and reduced oxygen extraction fraction). In an analogous fashion, the increase in anticorrelated DFA signal (decreasing Hurst exponent) seen in the 'blue' voxels should map decreased venous oxygenation or an increased oxygen extraction fraction. Follow-up studies relating 'blue' DFA maps to quantitative BOLD imaging or PET imaging with 15 O may permit more definitive interpretation. [26] Attenuation of the 'blue' DFA signal over time in both the square wave and ramp CO 2 sequences requires comment. A partial explanation may be surmised from the two studies done in the subjects with measures of cerebral oxygen saturation using near infrared spectroscopy and middle cerebral blood flow velocity with transcranial Doppler (Figures 6 and 7). In both these studies, cerebral blood flow velocity (as a surrogate for cerebral blood flow) showed a progressive increase in flow velocity that outlasted the maximal CO 2 response of the ramp sequence. In the first subject ( Figure 6) there was also a vigorous increase in blood pressure to the CO 2 stimulus. In the presence of increased CO 2 the cerebral circulation becomes pressure passive and cerebral blood flow would passively increase with increased driving pressure. The dynamic responsiveness, with resolution of 'blue' DFA voxels over time may be, in part, contingent on the systemic pressure mediated changes in response to CO 2 . Thus, to the extent that risk is related to either the regional changes in blood flow or its time course, the interplay between CVR and DFA responsiveness provides two windows from which to view individual patient neurological risk in the presence of steno-occlusive disease of the cerebral vasculature.
With DFA the histogram density peaked at a Hurst exponent near 1.0 with a decrease in the coefficient of variation. Healthy physiologic signals often are seen with a Hurst exponent approaching unity. Highlighting the one patient with both a square wave and ramp sequence suggests an adaptive response is seen to the perturbation of CO 2 stimulus. One interpretation from this patient is that these results support the concept of allometric control in response to a stimulus -measuring the fractal time signature of the signal response and its change over time. [23] West and colleagues have also advanced the notion of maximal information transfer in fractal physiological systems when signal noise is tuned to a fractal exponent of unity. [27] [28] In summary, we present two protocols to induce a CO 2 stimulus and two analysis techniques to process the resultant MR-BOLD signal. Such an approach demonstrates the potential utility of an enhanced brain stress test in the evaluation of patients with cerebrovascular steno-occlusive disease. Further studies based on these findings may delineate which patients are at greatest risk of stroke or post-operative delirium -to name but two examplesfollowing surgical procedures.

Supporting Information
Movie S1 Changes in CVR over time with the ramp stimulus in a patient with 'blue' brain. This video shows the alterations in CVR slope with ramp alterations in end-tidal CO 2 tensions. The dynamic behavior is in evidence here over the course of the provocative stimulus. (MP4)