# TUTORIAL 2: CVR mapping without respiratory challenges: i.e. using breath-holding or just a simple resting-state fMRI scan

This tutorial comprises a template data analysis script that utilizes native matlab functions as well as functions from the seeVR toolbox. If you use any part of this process or functions from this toolbox in your work please cite the following article:

Bhogal AA.,Medullary vein architecture modulates the white matter BOLD cerebrovascular reactivity signal response to CO2: observations from high-resolution T2* weighted imaging at 7T. Neuroimage 15 Dec, 2021. Vol 245, https://doi.org/10.1016/j.neuroimage.2021.118771

... and toolbox

Bhogal AA, (2021), abhogal-lab/seeVR: v1.53 (v1.53). Zenodo. https://doi.org/10.5281/zenodo.5283594

If you have suggestions for improvement, please contact me at a.bhogal@umcutrecht.nl, https://www.seevr.nl/contact-collaboration/ or feel free to start a new branch at https://github.com/abhogal-lab/seeVR

************************************************************************************************************************************

loadTimeseries: wrapper function to load nifti timeseries data

loadMask: wrapper function to load associated masks

meanTimeseries: calculates the average time-series signal intensity in a specific ROI

demeanData: subtracts the mean value from each voxel time-series

normTimeseries: normalizes time-series data to a specified baseline period

remLV: generates a mask that can be used to isolate and remove large vessel signal contributions

smthData/filterData: perfomes gaussian smoothing operation on image/time-series data

glmCVRidx: uses a general linear model approach to calculate a 'CVR index' - suitable when CO2 data is not available

lagCVR: calculates CVR and hemodynamic lag using a cross-correlation or lagged-GLM approach

### *************************************************************************************************************

### MRI Data Properties

This example uses BOLD data acquired using a Philips 7 tesla MRI scanner with the

following parameters:

Scan resolution (x, y): 112 109

Scan mode: MS

Repetition time [ms]: 1700

Echo time [ms]: 25

FOV (ap,fh,rl) [mm]: 224 59.2 224

Scan Duration [sec]: 281

Slices: 27

Volumes: 160

EPI factor: 37

Slice thickness: 2mm

Slice gap: 0.2mm

In-place resolution (x, y): 2mm

### MRI Data Processing

Simple pre-processing of MRI data was done using FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) as follows

1) motion correction (MCFLIRT)

2) calculate mean image (FSLMATHS -Tmean)

3) brain extraction (BET f = 0.2 -m)

4) tissue segmentation on BOLD image (FAST -t 2 -n 3 -H 0.2 -I 6 -l 20.0 -g --nopve)

# Analysis Pipeline

### 1) Setup options structure and load data

The seeVR functions share parameters using the 'opts' structure. In the current version this is implemented as a global variable (subect to change in the future). To get started, first initialize the global opts struct. Nifti data can be loaded using the loadTimeseries/loadMask wrapper functions.

clear all

clear global all

close all

% initialize the opts structure (!)

global opts

opts.verbose = 0;

% Set the location for the MRI data

datadir = 'ADDPATH\';

cd(datadir)

% Load motion corrected data

filename = 'BOLD_masked_mcf.nii.gz';

[BOLD,INFO] = loadTimeseries(datadir, filename);

% Load GM segmentation

file = 'BOLD_mean_brain_seg_0.nii.gz';

[GMmask,INFOmask] = loadMask(datadir, file);

% Load WM segmentation

file = 'BOLD_mean_brain_seg_1.nii.gz';

[WMmask,~] = loadMask(datadir, file);

% Load CSF segmentation

file = 'BOLD_mean_brain_seg_2.nii.gz';

[CSFmask,~] = loadMask(datadir, file);

% Load brain mask

file = 'BOLD_mean_brain_mask.nii.gz';

[WBmask,~] = loadMask(datadir, file);

### 2) Setup directories for saving output

Pay special attention to the savedir, resultsdir and figdir as you can use these to organize the various script outputs.

% specify the root directory to save results

opts.savedir = datadir;

% specify a sub directory to save parameter maps - this

% directory can be changed for multiple runs etc.

opts.resultsdir = fullfile(datadir,'RESULTS'); mkdir(opts.resultsdir);

% specify a sub directory to save certain figures - this

% directory can be changed for multiple runs etc.

opts.figdir = fullfile(datadir,'FIGURES'); mkdir(opts.figdir);

### 3) Take a quick look at the data

Use the 'meanTimeseries' function with any mask to look at the ROI timeseries. To see the differences between GM, WM and whole-brain timeseries, use the 'demeanData' function before plotting. We can also use the 'normTimesieres' function to normalize the BOLD data and look at the %BOLD signal change relative to the baseline period. There are two ways to do this: 1) specify the baseline start and end indexes; 2) select the baseline start and end points.

figure(10);

TS1 = meanTimeseries(BOLD, GMmask);

TS1 = demeanData(TS1);

% Plot GM timeseries

plot(TS1, 'k'); title('3X inspiratory breath-hold');

xlabel('absolute signal'); ylabel('time (s)'); xlim([0 length(TS1)]);

% Overlay WM timeseries

hold on

TS2 = meanTimeseries(BOLD, WMmask);

TS2 = demeanData(TS2);

plot(TS2, 'c');

title('3X inspiratory breath-hold'); xlim([0 length(TS1)]);

% Overlay WB timeseries

TS3 = meanTimeseries(BOLD, WBmask);

TS3 = demeanData(TS3);

plot(TS3, 'm');

title('3X inspiratory breath-hold'); xlim([0 length(TS1)]);

legend('GM','WM', 'WB')

hold off

clear TS1 TS2 TS3

### 4) Removing contributions from large vessels

CVR data can often be weighted by contributions from large vessels that may overshadow signals of interest. Particulary, for very high resolution acquisitions at high field strength. We can modify the whole brain mask to exclude CSF and large vessel contributions using the 'remLV' function ('remove large vessels'). This can be done by specifying the percentile threshold from which to remove voxels, or if the appropriate toolbox is not available, a manual threshold value.

% Supply necessary options

% Define the cutoff percentile (higher values are removed; default = 98);

opts.LVpercentile = 99;

% If the stats toolbox is not present, then a manual threshold must

% be applied. This can vary depending on the data (i.e. trial and error).

opts.verbose = 1;

[mWBmask] = remLV(BOLD,WBmask,opts);

saving tSD map
saving tSNR map
saving 1/tSD map
saving tNSR map
updating whole brain mask

### 5) Estimating relative CVR

In a perfect world we would have some idea of arterial gas tensions to which we could relate our BOLD signal changes. This would allow us to normalized our data and express signal changes as a function of changes in mmHg of CO2 or O2. For simple BH experiments this is often not the case - due to the need for gas monitoring systems etc. In cases where gas traces are not available, relative CVR maps can be generated using the 'glmCVRidx' function. This functions uses a reference signal (i.e. from a GM, cerebelum or sinus mask) as a explanatory variable and outputs a parameter map of beta values that can be interpreted as the relative CVR response. Signals are first band-pass filtered using the frequency range specified in the opts.fpass parameter. The function returns the CVRidx map, the band passed reference signal used for normalization and the band-passed BOLD timeseries if specified. Resulting maps are saved in opts.resultsdir\CVRidx\. NOTE this function can also be used to estimate CVR from resting state fMRI. In this case the vascular response is estimated based on low frequency oscillations. If using resting state data, it might be informative to explore different frequency limits for opts.fpass.

Often, breath-hold data can be somewhat corrupted due to motion artefacts. Imagine a subect is instructed to do an inspiratory BH (considered more tolerable than expiratory BH), when they inhale, the lungs expand and this can lead to tilting of the head. This motion can lead to artefacts in our BOLD signal that correlate with our BH task. We can examine the motion parameters and see if it is worthwhile to regress them out from our BOLD data. One of the strengths of the seeVR functions is the ability to play with explanatory regressors and nuissance regressors. In tutorial 1 this is done using the scrubData.m function. In the current seeVR version, handling of nuissance signals has been integrated in the glmCVRidx function.

% Normalize to baseline using interactive mode

baseline_index = [1 20];

nTS1 = normTimeseries(BOLD, GMmask, baseline_index);

opts.spatialdim = 3; % specify 3D gaussian smoothing (default 2D smoothing)

opts.FWHM = [3 3 3]; % specify FWHM of smoothing kernel (default 4mm)

% 'standard' gaussian smoothing - also see filterData function below for more

% smoothing options

sBOLD = smthData(nTS1, WBmask, opts);

Finished smoothing Data

clear nTS1

% TRY EDGE PRESERVING SMOOTHING! guided (edge preserving) smoothing using bilateral filter

%opts.filter = 'gaussian' or 'tips' or 'bilateral

%opts.sigma_range = 20;

% this is a tunable parameter for the bilater filter. Higher values prodice a more 'gaussian'-like kernel.

% Lower values do better at preserving edges (i.e. between GM/WM

% boundaries). I have no 'right' way to apply this; data-dependent + trial

% & error

% opts.FWHM = 8; % specify FWHM of smoothing kernel (default 4mm). This

% usually needs to be larger in order to include enough voxels

%guideim = squeeze(BOLD(:,:,:,1)); %use first BOLD image as quide for smoothing

%[sBOLD] = filterData(nBOLD,guideim,mWBmask,opts);

% to calculate CVR, nuissance regressors can now be supplied.

% ** note, this function can also output band-passed signals

% ** note, here we supply the GMmask as an ROI to generate the reference

% signal but you can supply any mask (for example, cerebellum or CSF or a

% saggital sinus mask etc.)

% load nuisance signals

% load nuisance signals

cd(datadir) % Go to our data directory

mpfilename = 'BOLD_mcf.par'; % Find MCFLIRT motion parameter file

nuisance = load(mpfilename); % Load nuisance regressors (translation, rotation etc.)

% Calculate motion derivatives

dtnuisance = gradient(nuisance); % temporal derivative

sqnuisance = nuisance.*nuisance; % square of motion

motion = [nuisance dtnuisance sqnuisance];

drift_term = 1:1:size(sBOLD,4);

%concatenate motion params with drift term

np = [motion drift_term'];

% Demean/rescale regressors

for ii=1:(size(np,2)); np(:,ii) = demeanData(rescale(np(:,ii))); end

opts.fpass = [0.000001 0.15]; % accept a wide frequency band (for LFO you can set this to 0.0001 --> 0.08 for example

[CVRidx, ~, ~] = glmCVRidx(sBOLD, WBmask, GMmask, np, opts);

% and now without nuisance regression

%[CVRidx2, ~,~] = glmCVRidx(sBOLD, WBmask, GMmask, [], opts);

% plot the effect of nuisance signal regression on CVR

figure;

step = 5;

z=4;

scl = [0 7];

for ii=0:step-1

% with

subplot(1,5,ii+1)

imagesc(imrotate(CVRidx.CVRidx(:,:,z + ii*step),90), scl);

colormap(flip(brewermap(256,'Spectral')));

freezeColors;

cleanPlot; if ii==0; ylabel('CVR with nuisance', 'color', 'w'); end

end

set(gcf, 'Units', 'pixels', 'Position', [10, 80, 1000, 200]);

clear nTS1 %clear large matrices to free up data

### This approach outlined above can also be used to generate CVR maps from resting-state fMRI data!

### 6) Calculating Hemodynamic Lag

In addition to CVR, hemodynamic lag can be another interesting parameter to consider when looking at vascular impairment. A common approach to doing this involves performing a cross-correlation or lagged-GLM approach with some stimulus waveform (PetCO2 or PetO2). The calculated lags with respect to this waveform probe are then plotted voxel-wise to produce a lag parameter map. A more advanced approach is to use this input probe in order to generate an optimized regressor (see: https://rapidtide.readthedocs.io/en/stable/). This analysis technique has also been implemented in the seeVR toolbox in the 'lagCVR' function. This function is capable of generating an optimized BOLD regressor, lag maps based on both a cross-correlation and/or lagged-GLM analysis, CVR maps as well as many statistical maps to evaluate analysis outputs. Since we do not have any breathing traces in this example, we will use the signal as a probe through which to select highly correlated voxels that will comprise the optimized BOLD regressor. Furthermore, we will skip generating CVR maps (lagCVR performs linear regression of the BOLD signal against the end-tidal trace to generate CVR/O2-resonse maps) - we've anyways already created the CVR index maps in step 5. As input, use the band-passed signals obtained using glmCVRidx above.

% First lets smooth our cleaned BOLD (step 4) data using the same smoothing

% parameters applied in step 6. We're using a lot of smoothing mostly for

% visual purposes - the SNR of 7T BOLD data is generally quite high. Use

% the modified WBmask from step 5 to avoid some of the larger vessels.

% The lagCVR function has many options and possible outputs; some of which

% more useful than others. For details, look at default options in the .m file.

% lagCVR has the option to add a nuissance regressor (for example a PetCO2

% trace when focusing on using PetO2 as the main explanatory regressor).

% This regressor is only used when opts.glm_model = 1 (default = 0). For

% this example lets run the GLM and include the original motion parameters

% and drift term generated above. Since we focus on the GLM, lets skip the correlation

% approach all together.

opts.glm_model = 1;

opts.corr_model = 1;

% an alternate way to produce CVR maps - these should be similar to those

% derived in the step above

opts.cvr_maps = 1; %default = 1

% Produces the optimized regressor (default = 1). When this is set to 0; a straight correlation

% with the input probe is done

opts.refine_regressor = 1; %default = 1

% Factor by which to temporally interpolate data. Better for picking up

% lags between TR. Higher value means longer processing time and more RAM

% used (so be careful) - set to 1 for fastest results!

opts.interp_factor = 2; %default is 4

% The correlation threshold is an important parameter for generating the

% optimized regressor. For noisy data, a too low value will throw an error.

% Ideally this should be set as high as possible, but may need some trial

% and error.

opts.corr_thresh = 0.6; %default is 0.7

% Thresholds for refinening optimized regressors. If you make this range too large

% it smears out your regressor and lag resolution is lost. When using a CO2

% probe, the initial 'bulk' alignment becomes important here as well. A bad

% alignment will mean this range is also maybe not appropriate or should be

% widened for best results. Since we use the actual BOLD as a probe we keep

% this range tight. (ASSUMING TR ~1s!, careful).

opts.lowerlagthresh = -2; %default is -3

opts.upperlagthresh = 2; %default is 3

% Lag thresholds for lag map creation. Since we are looking at a healthy

% brain, this can be limited to between 20-60TRs. For impariment you can consider

% to raise the opper threshold to between 60-90TRs. (ASSUMING TR ~1s!, careful).

opts.lowlag = -5; %setup lower lag limit; negative for misalignment and noisy correlation

opts.highlag = 30; %setups upper lag limit; allow for long lags associated with pathology

Perform hemodyamic analysis

% For the GLM analysis, it is possible to supply nuisance parameters. Try to

% compare results using the spatially smoothed (sBOLD) data with and

% without nuisance signals (np) - BEWARE, results are not always better...

probe = meanTimeseries(sBOLD, GMmask); % reference signal

[optimized_probe, maps] = lagCVR(GMmask,mWBmask,sBOLD,probe,drift_term,opts);

probename = 'final_probe.mat'

Creating optimized regressor
Correlating TS voxelwise to probe 1
.
.
.
Elapsed time is 0.046261 seconds.
finished correlation, estimating new probe

rmse = 0.2653

Correlating TS voxelwise to probe 2
.
.
.
Elapsed time is 0.040881 seconds.
finished correlation, estimating new probe

rmse = 0.0035

Finished creating optimized regressor
performing correlation based lag analysis using optimized regressor

passes = 1

perc = 0.0157

0 percent of voxels have clipped lag values
performing correlation based lag analysis using input probe

passes = 1

perc = 0.0470

0 percent of voxels have clipped lag values
Lag, regressor and r_maps were created in: 0 minutes
Generating base maps using probe trace
Base CVR, r2, Tstat and SSE maps were created using entidal regressor
Generating lag-adjusted maps
Stimulus response data is saved
performing GLM based lag analysis using OPTIMIZED regressor
Generating lag-adjusted CVR maps based on GLM analysis
performing GLM based lag analysis using INPUT regressor
Generating lag-adjusted CVR maps based on GLM analysis
finished running GLM analysis in: 2 minutes
saving maps in .mat file

[~, maps2] = lagCVR(GMmask,mWBmask,sBOLD,probe,np,opts);

probename = 'final_probe.mat'

Creating optimized regressor
Correlating TS voxelwise to probe 1
.
.
.
Elapsed time is 0.027345 seconds.
finished correlation, estimating new probe

rmse = 0.2653

Correlating TS voxelwise to probe 2
.
.
.
Elapsed time is 0.037731 seconds.
finished correlation, estimating new probe

rmse = 0.0035

Finished creating optimized regressor
performing correlation based lag analysis using optimized regressor

passes = 1

perc = 0.0157

0 percent of voxels have clipped lag values
performing correlation based lag analysis using input probe

passes = 1

perc = 0.0470

0 percent of voxels have clipped lag values
Lag, regressor and r_maps were created in: 0 minutes
Generating base maps using probe trace
Base CVR, r2, Tstat and SSE maps were created using entidal regressor
Generating lag-adjusted maps
Stimulus response data is saved
performing GLM based lag analysis using OPTIMIZED regressor
Generating lag-adjusted CVR maps based on GLM analysis
performing GLM based lag analysis using INPUT regressor
Generating lag-adjusted CVR maps based on GLM analysis
finished running GLM analysis in: 2 minutes
saving maps in .mat file

% The advantage of using the global opts struct is that the variables used

% for a particular processing run (including all defaults set within

% functions themselves) can be saved to compare between runs.

save([opts.resultsdir,'processing_options.mat'], 'opts');

### 7) Plot Results

In the maps struct, we have a subfield called 'XCORR' where the results of the cross-correlation based lag results are returned. Along with using the optimized BOLD regressor, in this case we also performed a straight correlation with the input proble (this can be turned off using opts.trace_corr = 0). When both the optimized regressor and input regressor are used, the maximum r value compared between each analysis approach is used to define the lag. This is useful if one regressor is more accurate than the other for a particular voxel. Generally though, the optimized regressor (lag_opti) provides the best results.

% Plot resulting lag map with CVR maps

lagMap1 = maps.GLM.optiReg_lags;

lagMap2 = maps2.GLM.optiReg_lags;

lagMap1(lagMap1 == 0) = NaN;

lagMap2(lagMap2 == 0) = NaN;

CVR = CVRidx.CVRidx;

figure;

step = 5;

z=2;

sclCVR = [0 7];

scllag = [3 15];

for ii=0:step-1

% BOLD

subplot(4,5,ii+1)

imagesc(imrotate(BOLD(:,:,z + ii*step,1),90));

colormap(gray)

freezeColors;

cleanPlot; if ii == 0; ylabel('BOLD', 'color', 'w'); end

%CVR

subplot(4,5,(ii+1+step))

imagesc(imrotate(CVR(:,:,z + ii*step),90), sclCVR);

colormap(viridis);

freezeColors;

cleanPlot; if ii == 0; ylabel('CVR', 'color', 'w');

xlabel('relative CVR: 0 to 10', 'color', 'w'); end

%lags

subplot(4,5,(ii+1+2*step))

imagesc(imrotate(lagMap1(:,:,z + ii*step),90), scllag);

%colormap(flip(brewermap(256,'Spectral')))

colormap(magma)

freezeColors;

cleanPlot; if ii == 0; ylabel('Lags without nuisance', 'color', 'w');

xlabel('lags: 3 to 15s', 'color', 'w'); end

subplot(4,5,(ii+1+3*step))

imagesc(imrotate(lagMap2(:,:,z + ii*step),90), scllag);

%colormap(devon)

colormap(magma)

freezeColors;

cleanPlot; if ii == 0; ylabel('Lags with nuisance', 'color', 'w');

xlabel('lags: 3 to 15s', 'color', 'w'); end

end

set(gcf, 'Units', 'pixels', 'Position', [10, 100, 1000, 840]);