The test data for this tutorial consists of a pre-processed BOLD dataset acquired using a 7 tesla Philips MRI system with corresponding end-tidal CO2 traces measured using a 3rd generation RespirAct computer controlled gas delivery system (https://thornhillmedical.com/research/respiract-ra-mr/). The respiratory protocol consisted of a baseline period followed by a hypercapnic block and then a final baseline period. Unzip the data file and update the data path (/Data) in the included Matlab live script (.mlx) file.

TUTORIAL 3: Dispersion mapping as an alternative to hemodynamic lag

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

Alex A. Bhogal, (2023). abhogal-lab/seeVR: 1.53 (v1.53). Zenodo. https://doi.org/10.5281/zenodo.7816690

If you have suggestions for improvement or find bugs, 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 using native matlab functions

loadMask: wrapper function to load nifti mask/image data using native matlab functions

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

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

denoiseData: temporally de-noises data using a wavelet or moving-window based method (toolbox dependent)

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

lagCVR: calculates CVR and hemodynamic lag using a cross-correlation or lagged-GLM approach. Hemodynamic maps and optimized BOLD regressor are output.

fitTau: fits the dispersion time constant and associated model parameters

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

MRI Data Properties

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

following parameters:

Scan resolution (x, y): 136 133

Scan mode: MS

Repetition time [ms]: 3000

Echo time [ms]: 25

FOV (ap,fh,rl) [mm]: 217.600 68.800 192.000

Scan Duration [sec]: 369

Slices: 43

Volumes: 120

EPI factor: 47

Slice thickness: 1.6mm

Slice gap: none

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

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 -o)

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

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 based on the native matlab nifti functions (or 'loadImageData' using provided nifti functions - see which one works for you). This function will also initialize certain variables in the opts structure including opts.TR (repetition time). opts.dyn (number of volumes), opts.voxelsize (resolution), and opts.info ( if using loadTimeseries then opts.headers is initialized and used by saveImageData to save timeseries, parameter maps and masks in opts.headers.ts/map/mask).

*If you use your own functions to load your imaging data, ensure that you also fill the above fields in the opts structure to ensure maintained functionality (especially opts.TR). If functions throw errors, it is usually because a necessary option is not specified before the function call.

Startup - use 'ctrl + enter' to run individual code blocks

 
clear all
clear global all
close all
 
% initialize the opts structure (!)
global opts
addpath(genpath('ADDPATH TO seeVR toolbox'));
 
% Set the location for the MRI data
datadir = 'ADDPATH to DATA';
 
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.

Setup necessary directories

 
% specify the root directory to save results
opts.savedir = datadir; 
 
% specify a results 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 'normTimeseries' function before plotting using the 'meanTimeseries' function. As you can see in the figure, there are three signal peaks corresponding to three hypercapnic periods delivered using a RespirAct system.

Visualize ROI timeseries

% Normalized BOLD data to 20 volumes in baseline period

% if no index is provided, baseline can be selected manually

nBOLD = normTimeseries(BOLD, WBmask, [5 25]);

TS1 = meanTimeseries(nBOLD, GMmask);

%establish new xdata vector for plotting

xdata = opts.TR:opts.TR:opts.TR*size(nBOLD,4);

% Plot GM timeseries

figure;

plot(xdata, TS1, 'k'); title('time-series data');

ylabel('absolute signal'); xlabel('time (s)'); xlim([0 xdata(end)])

% Overlay WM timeseries

hold on

TS2 = meanTimeseries(nBOLD, WMmask);

plot(xdata, TS2, 'm');

set(gcf, 'Units', 'pixels', 'Position', [200, 500, 600, 160]);

legend('GM','WM')

hold off

clear TS1 TS2 nBOLD

4: Load end-tidal gas traces

For simplicity I have already aligned the gas traces. Simply load them from the example dataset. If loading your own data, you can use the 'resampletoTR' function to interpolate breathing data to the same temporal resolution as your MRI data. For alignment, use the trAlign function (see tutorial 1). If you need help getting your physiological data loaded in, email me and I can write you a function.

cd(datadir)
load('traces.mat')

5: Removing contributions using a large vessel mask

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 = 98;

% 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).

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

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

6: Temporally de-noise data

To remove spikes from motion or high frequency noise these try the 'denoiseData' function. When the wavelet toolbox is available, this function uses a wavelet-based approach. Otherwise a more simple moving window approach is applied, however this can result in a temporal shift. Low-pass filtering is another option - this can be done using the bandpassfilt.m function.

Temporal de-noising

%If no wavelet toolbox, then moving average is applied. See

%function for options.

opts.wdlevel = 2; %level of de-noising (higher number = greater effect)

opts.family = 'db4'; %family - can opitimize based on data

opts.verbose = 1;

denBOLD = denoiseData(BOLD, WBmask, opts);

finished discreet wavelet denoising in: 6 seconds

7: Smooth and normalize data

There are several tunable smoothing options ranging from standard gaussian to edge-preserving bilateral smoothing. For MAC users (sorry) this is restricted to gaussian but you can easily replace with your own spatial smoothing algorithms.

Spatial smoothing

 
% 'guided' gaussian smoothing (for MAC see smthData function)
opts.filter = 'gaussian';
opts.spatialdim = 3; %default = 2
opts.FWHM = [5 5 5];
 
% Normalize data to first 15 baseline images. *NB if nornIdx is not
% supplied, you will be asked to manually select baseline indices.
normIdx = [1 15];
nBOLD = normTimeseries(denBOLD,mWBmask,normIdx);
 
%use modified mask from step 5 as input to avoid smoothing large vessel
%signals into tissues
guideImg = squeeze(BOLD(:,:,:,1));
sBOLD = filterData(nBOLD, guideImg, mWBmask, opts);
 
clear nBOLD denBOLD

8) Dispersion Mapping (including lag mapping for comparison)

Hemodynamic lag calculations (see tutorial #1) provides information on possible blood flow delays. Generally, the lag is expressed in seconds or TR but is not always reflective of true temporal effects. Often, correlations can be weighed by the time-to-peak of the BOLD signal response. This can be seen when comparing white matter and grey matter time courses. The white matter signal rises more slowly and this leads to longer delays. However, the onset of the white matter response is often close to that of the GM. This effect is can be termed 'signal dispersion', and may be related to draining vein effects. See: https://doi.org/10.1016/j.neuroimage.2021.118771 and https://pubmed.ncbi.nlm.nih.gov/26126862/.

The 'fitTau.m' function can be used to calculate the voxel-wise dispersion. This function also outputs a scaling map that can be considered as a 'disperson-corrected' CVR map.

%First we will generate an initial probe using the lagCVR function. 
%This probe will represent the 'fastest' intrinsic responses to the 
%vasoactive stimulus and will already show inherent changes that 
% occur as the CO2 bolus moves from the lungs to the brain and then 
 
opts.corr_model = 1; 
opts.glm_model = 1; %this can take a while - to explore and increase processing speed just use the corr_model. NB that turning this off will cause problems with some plots below.
opts.cvr_maps = 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)
opts.interp_factor = 1; %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.corrthresh = 0.7; %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. (ASSUMING TR ~1s!, careful). 
opts.lowerlagthresh = -2; %default is -3
opts.upperlagthresh = 2; %default is 3
 
% Lag thresholds (in units of TR) 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 = -3; %setup lower lag limit; negative for misalignment and noisy correlation
opts.highlag = 25; %setups upper lag limit; allow for long lags associated with pathology
 
% Perform hemodyamic analysis
% The lagCVR function saves all maps and also returns them in a struct for
% further analysis. It also returns the optimized probe when applicable.
 

Lets load our motion parameters to compare the effect of adding them for lag mapping

cd(datadir) % Go to our data directory

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

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

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

%concatenate motion params with drift term

np = [nuisance drift_term'];

[newprobe, maps] = lagCVR(GMmask, mWBmask, sBOLD, CO2trace, np, opts);

using input nuisance signals and temporal derivative
adding drift term
rescaling nuisance between -1 & 1
filtering regressors based on correlation value supplied in opts.motioncorr

probename = 'final_probe.mat'

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

rmse = 1.2101e+03

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

rmse = 0.0092

Correlating TS voxelwise to probe 3
 . 
 . 
 . 
finished correlation, estimating new probe
Elapsed time is 2.364854 seconds.

rmse = 0.0057

Correlating TS voxelwise to probe 4
 . 
 . 
 . 
finished correlation, estimating new probe
Elapsed time is 2.446702 seconds.

rmse = 0.0041

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

passes = 1

perc = 4.9737

5 percent of voxels have clipped lag values

passes = 2

perc = 2.6942

3 percent of voxels have clipped lag values

passes = 3

perc = 2.0019

2 percent of voxels have clipped lag values

passes = 4

perc = 1.6810

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

passes = 1

perc = 8.8149

9 percent of voxels have clipped lag values

passes = 2

perc = 4.7155

5 percent of voxels have clipped lag values

passes = 3

perc = 3.7246

4 percent of voxels have clipped lag values

passes = 4

perc = 3.3259

3 percent of voxels have clipped lag values

passes = 5

perc = 3.1388

3 percent of voxels have clipped lag values

passes = 6

perc = 3.0378

3 percent of voxels have clipped lag values

passes = 7

perc = 2.9726

3 percent of voxels have clipped lag values

passes = 8

perc = 2.9285

3 percent of voxels have clipped lag values

passes = 9

perc = 2.9049

3 percent of voxels have clipped lag values

passes = 10

perc = 2.8881

3 percent of voxels have clipped lag values
exceeded the maximum allowed passes
... to increase passes set opts.passes to a higher value
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 endtidal 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

Fitting Dispersion

Dispersion fitting has been simplified compared the the older tutorial. The option to generate a look-up table including onset is still possible using the convHRF2 and firHRF2 functions.

 
% 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.
 
[tau_maps] = fitTau(CO2trace, sBOLD, mWBmask, opts)
passes = 1
3 percent of voxels have clipped tau values
passes = 2
2 percent of voxels have clipped tau values
passes = 3
2 percent of voxels have clipped tau values
passes = 4
2 percent of voxels have clipped tau values
Elapsed time is 150.435448 seconds.
tau_maps = struct with fields:
    expHRF: [1×1 struct]
 
save([opts.resultsdir,'processing_options.mat'], 'opts');

9: Plot results

Compare CVR with Lag-Corrected CVR

%setup options

opts.scale = [-0.4 0.4]; % this is the expected data range. Default is [-5 5]

opts.row = 5; % this is the nr of rows. More rows means more images. Too many images will throw an error

opts.col = 6; % this is the nr of columns. This should be an even number. For 6 cols, 3 will be source images and 3 will be the param map.

opts.step = 1; % This value is multiplied by 2 in the function. So step = 2 means you jump 4 slices between images.

opts.start = 8; % this is the starting image

%PLOT CVR

%rotate input images to display correctly

sourceImg = imrotate(guideImg,90); %use the guide image we used for smoothing. Alternatively a mean BOLD image or single time-point

paramMap1 = imrotate(maps.XCORR.CVR.bCVR, 90); % use the basic CVR map calculated using LagCVR

mask = imrotate(mWBmask, 90);

map = flip(brewermap(128, 'Spectral')); %use the spectral colormap (or any other you like)

%generate plots

plotMap(sourceImg,mask,paramMap1,map,opts);

%PLOT LAG-CORRECTED CVR

paramMap2 = imrotate(maps.XCORR.CVR.cCVR, 90);

plotMap(sourceImg,mask,paramMap2,map,opts);

%PLOT DIFFERENCE - high values indicate where CVR estimate may have been

%improved after considering lag

opts.scale = [-0.02 0.02];

map = BWR2(51);

plotMap(sourceImg,mask,(paramMap2 - paramMap1),map,opts);

Compare LAG with TAU map

%PLOT HEMODYNAMIC LAG MAP USING OPTIMIZED BOLD REGRESSOR AND NUISANCE

%SIGNALS (GLM MODEL)

opts.scale = [0 60];

%update colormap

map = jet(60);

paramMap1 = imrotate(maps.GLM.optiReg_lags, 90); % use the basic CVR map calculated using LagCVR

%plot map