TUTORIAL 1: BOLD-CVR based on computer controlled CO2 stimulus
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:
... and toolbox
***********************************************************************************************************************
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
loadRAMRgen4: loads respiratory traces generated using a 4th generation RespirAct system
trAlign: aligns respiratory traces with some input signal
- alternate function 'autoAlign' that works most optimally using fixed inspired systems (less sharp transitions = less error)
- combine with 'resampletoTR' when loading physiological trace text files
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)
smthData: performs gaussian smoothing operation on image/time-series data (also edge-preserving options - see function)
- alternate function 'filterData' for different smoothing options (not implemented for MAC)
lagCVR: calculates CVR and hemodynamic lag using a cross-correlation or lagged-GLM approach
*************************************************************************************************************
MRI Data Properties
The protocol applied here consisted a longer hypercapnic period followed by two shorter hypercapnic blocks. BOLD data was acquired using a Philips 3 T scanner during manipulation of arterial CO2 and O2 gases:
Scan resolution (x, y): 88 88
Scan mode: MS (MULTI-BAND)
Repetition time [ms]: 1050
Echo time [ms]: 30
FOV (ap,fh,rl) [mm]: 220 128 220
Slices: 27
Volumes: 496
EPI factor: 49
Slice thickness: 2.5mm
Slice gap: none
In-plane resolution (x, y): 2.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) perform distortion correction (TOPUP)
3) calculate mean image (FSLMATHS -Tmean)
4) brain extraction (BET f = 0.2 -m)
5) 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
% initialize the opts structure (!)
%add seeVR to Matlab path
addpath(genpath('ADDPATH TO seeVR toolbox'));
% Set the location for the MRI data
datadir = 'ADDPATH to DATA\';
% Set the location for the end-tidal data
opts.seqpath = 'ADDPATH to DATA\traces';
%!!!! if input data is type int16/int32 (sometimes in SIEMENS data) then
%use the loadImagedata function for all images including masks!!!!
% Load motion corrected data
filename = 'BOLD_applytopup.nii.gz';
[BOLD,INFO] = loadTimeseries(datadir, filename);
file = 'BOLD_mean_brain_seg_0.nii.gz';
[GMmask,INFOmask] = loadMask(datadir, file);
file = 'BOLD_mean_brain_seg_1.nii.gz';
[WMmask,~] = loadMask(datadir, file);
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
% 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(xdata, TS1, 'k'); title('time-series data');
ylabel('absolute signal'); xlabel('time (s)'); xlim([0 xdata(end)])
TS2 = meanTimeseries(nBOLD, WMmask);
set(gcf, 'Units', 'pixels', 'Position', [200, 500, 600, 160]);
4: Load end-tidal gas traces and temporally align with MRI data
For this purpose you can use scripts located in the seeVR/gas_import/ folder (these are 'work in progress'). Typically, cross correlation between the gas trace of interest and a whole brain or gray-matter signal is used for initial bulk alignment. This approach is prone to errors due to diffrences in the end expired gas trace vs. the hemodynamic response, as well as noise or other artefacts that can influence the correlation. For this reason, seeVR includes an alignment tool. If you are using your own loading scripts for the end-tidal traces ensure that they are first resampled to the TR of your MRI data. This can be done using the 'resampletoTR' function included in the gas_import folder. You can also try the autoAlign.m function which seems to work quite OK when using fixed inspired CO2 systems.
Align physiological traces
% Load breathing data based on RespirAct system
[CO2_corr,O2_corr,Respiration_rate_corr] = loadRAMRgen4(opts);
% Use the GM time-series for alignment
TS = meanTimeseries(BOLD,GMmask);
% Run GUI to correct alignment. If you dont have O2 data, you can supply
% only the CO2 data. Correlation with the time-series is done using the first input.
trAlign(CO2_corr, O2_corr, TS, opts); % for this example the alignment offset was 8 which is saved in the workspace
% If both CO2 & O2 are supplied, the GUI returns two probe vectors to the
% Matlab workspace. If only CO2 or O2 are provided, then one probe is
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 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. If this function fails, try restarting fromt step 1 and loading all data using loadImageData.m function.
Remove large vessel contributions
% Supply necessary options
% Define the cutoff percentile (higher values are removed; default = 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
opts.wdlevel = 2; %level of de-noising (higher number = greater effect)
opts.family = 'db4'; %family - can opitimize based on data
denBOLD = denoiseData(BOLD, WBmask, opts);
finished discreet wavelet denoising in: 9 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
% Normalize data to first 15 baseline images. *NB if nornIdx is not
% supplied, you will be asked to manually select baseline indices.
nBOLD = normTimeseries(denBOLD,mWBmask,normIdx);
%use modified mask from step 5 as input to avoid smoothing large vessel
guideImg = squeeze(BOLD(:,:,:,1));
sBOLD = filterData(nBOLD, guideImg, mWBmask, opts);
8: Calculate CVR and hemodynamic lag
Finally on to the fun stuff. The lagCVR function takes the PetCO2 input (or PetO2) and returns CVR maps and hemodynamic lag maps along with several statistical maps for data evaluation. This function has many options that can all impact your results. I've tried to set some basic defaults that should work with a variety of input data. For more information see the manual or explore the code and just try stuff.
Calculating CVR and Lag
% lagCVR has the option to add a nuissance regressor (for example a PetCO2
% trace when focusing on using PetO2 as the main explanatory regressor).
% Setup some extra options - most are default but we will initialized them
% for the sake of this tutorial
% We would like to generate CVR maps
opts.cvr_maps = 1; %default is 1
% Factor by which to temporally interpolate data. Better for picking up
% lags between TR. Higher value means longer processing time and more RAM
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
opts.corr_thresh = 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 = -5; %setup lower lag limit; negative for misalignment and noisy correlation
opts.highlag = 40; %setups upper lag limit; allow for long lags associated with pathology
% For comparison we can also run the lagged GLM analysis: beware this can
% take quite some time and memory
opts.glm_model = 1; %default is 0
opts.corr_model = 1; %default is 1
%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_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 3.606511 seconds.
rmse = 2.1149e+03
Correlating TS voxelwise to probe 2
.
.
.
finished correlation, estimating new probe
Elapsed time is 3.202075 seconds.
rmse = 0.0047
Finished creating optimized regressor
performing correlation based lag analysis using optimized regressor
passes = 1
perc = 9.9379
10 percent of voxels have clipped lag values
passes = 2
perc = 4.4259
4 percent of voxels have clipped lag values
passes = 3
perc = 2.5423
3 percent of voxels have clipped lag values
passes = 4
perc = 1.8660
2 percent of voxels have clipped lag values
performing correlation based lag analysis using input probe
passes = 1
perc = 10.2802
10 percent of voxels have clipped lag values
passes = 2
perc = 2.3431
2 percent of voxels have clipped lag values
passes = 3
perc = 1.4190
1 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 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: 0 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');
9: Plot results
Compare CVR with Lag-Corrected CVR
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 = 18; % this is the starting image
%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)
plotMap(sourceImg,mask,paramMap1,map,opts);
paramMap2 = imrotate(maps.XCORR.CVR.cCVR, 90);
plotMap(sourceImg,mask,paramMap2,map,opts);
opts.scale = [-0.02 0.02];
plotMap(sourceImg,mask,(paramMap2 - paramMap1),map,opts);
Compare LAG with and without nuisance
%PLOT HEMODYNAMIC LAG MAP USING OPTIMIZED BOLD REGRESSOR WITHOUT NUISANCE
paramMap1 = imrotate(maps.XCORR.lag_opti, 90); % use the basic CVR map calculated using LagCVR
plotMap(sourceImg,mask,paramMap1,map,opts);