2.5 Estimation of multiple hemodynamic maps from the concomitant CO2 and O2 data

Data analysis of the concomitant CO₂ and O₂ data was conducted using the software Statistical Parametric Mapping (SPM, University College London, UK) and in-house MATLAB (MathWorks, Natick, MA) scripts. In pre-processing, motion correction was performed by realigning the image volumes of the BOLD MRI images to the first volume within the time series. The BOLD images were then normalized to the image template of Montreal Neurological Institute (MNI) via T₁-weighted high-resolution anatomic image and were resampled to a voxel size of 2×2×2 mm³. As the final step of the pre-processing, all BOLD image volumes were smoothed using a Gaussian filter with a full-width-half-maximum (FWHM) of 8mm.

The post-processing followed the steps outlined in Figure 2. Step 1: estimation of bolus arrival time (BAT). BAT of the inhaled gas was quantified as the time delay between the end-tidal time courses (which has been corrected for sampling tubing delay) and the BOLD signal time courses. We have tested several approaches for this calculation and propose the following procedure for best performance. We used BOLD time course in a reference brain region as an intermediate and calculate voxel-by-voxel delay time between each brain voxels and the reference time course. For healthy volunteers, the reference region was chosen to be whole brain gray matter. For Moyamoya patients, since anterior and middle cerebral arteries could be affected but posterior cerebral circulation is thought to be intact in this disease, cerebellum gray matter time course was used as the reference. The best delay was determined by shifting the reference time course from −10·TR to +80·TR at an interval of TR and, at each shift step, performing linear regression between the voxel time course and the shifted reference. The shift that corresponds to the least residual signal in the regression analysis will be used as the best delay. Next, the delay between the reference time course and the end-tidal time course were calculated in a similar fashion. This was performed separately for EtCO₂ and EtO₂ time courses. The voxel-by-voxel BAT was then obtained by the sum of the voxel-vs-reference delay and the reference-vs-end-tidal time course delay. One therefore obtains a CO₂ BAT and an O₂ BAT map, although their image contrast is the same.

Illustration of analysis steps in the multiparametric imaging method to obtain bolus arrival time (BAT), cerebrovascular reactivity (CVR), venous cerebral blood volume (vCBV) and functional connectivity maps.

Step 2: estimation of CVR and vCBV (Figure 2). After obtaining the BAT maps, the EtCO₂ and EtO₂ time courses were shifted for each voxel. Then voxel-wise regression was performed, with the BOLD time course of each voxel as dependent variable and the BAT-shifted end-tidal time courses as independent variable. The magnitudes of CVR and vCBV of each voxel were calculated using the coefficients of the linear regression (Thomas et al., 2014; Yezhuvath et al., 2009), i.e., CVR(i, j, k) = b(i, j, k)/[c(i, j, k) + min(EtCO₂) · b(i, j, k)] and vCBV(i, j, k) = a(i, j, k)/[c(i, j, k) + min(EtO₂) · a(i, j, k)], and is written in %ΔBOLD signal per mmHg CO₂ and %ΔBOLD signal per mmHg O₂, respectively.

Step 3: estimation of functional connectivity maps (Figure 2). For functional connectivity analysis, the residual BOLD time course after factoring out the CO₂ and O₂ influences, i.e., ΔBOLD/BOLD₀ – a · EtO₂ – b · EtCO₂, was first calculated on a voxel-by-voxel basis, which forms a new 4D dataset. Then the residual 4D data were detrended and bandpass-filtered to 0.01–0.1 Hz to retain the low-frequency fluctuation components, which represent the main signal for functional connectivity. Next, independent component analysis (ICA) was performed using MELODIC (FMRIB Analysis Group, Oxford University). Group ICA analysis was used for the healthy volunteers, and individual ICA analysis was performed for Moyamoya patients. The component number of the ICA analysis was set to be 20, following previous literature (Lerman et al., 2014). The brain networks were identified by evaluating the similarity in the spatial extent between the ICA components and previously published results (Beckmann et al., 2005) using spatial cross-correlation. Visual inspection was performed to verify the identified networks.