Prep fMRI data for connectivity or computational model-based anlayses
Below is a step‐by‐step tutorial for pre‐processing fMRI data for connectivity analyses and custom linear or computational reinforcement learning (RL) models using CANLAB tools. In this example, we create a denoised 4-D time series (an fmri_data object) for each participant, then use it for downstream analyses. The CANLAB toolboxes (see CANLAB.github.io) offer functions for most of these steps. Feel free to adapt the code snippets to your specific dataset and analysis goals.
See also:
canlab_connectivity_preproc()
Table of Contents
A brief overview of steps
Step 1. Set Up Your Environment
Step 2. Create the fmri_data Object
Sample data
Orientation to Pinel data files
Step 3. Build the Nuisance Covariate Model
Step 4. Regress Out Nuisance Signals
Step 5. Save the Preprocessed Data
Step 6. Optional: Explore how nuisance covs overlap with task design
Step 6. Analysis: Connectivity and Custom RL Models
Summary
% Some formatting stuff for html
set(groot, 'DefaultFigurePosition', [200, 200, 900, 600], 'DefaultAxesFontSize', 12);
A brief overview of steps
The code below loads some sample data and does the complete filtering and nuisance regression
% Load the key 4-D image file into an object
fname = which('swrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz');
obj = fmri_data(fname);
% Obtain a matrix of 'spike' indicators and movement-related regressors
% based on multiple criteria to use as covariates
% 24 movement-related covariates per run
mvmtfname = which('rp_sub-sid001567_task-pinel_acq-s1p2_run-03_bold.txt');
[mvmt_matrix, mvmt_regs_24] = movement_regressors(mvmtfname);
covs = mvmt_regs_24;
% 'spike' outliers and high framewise-displacement (fd) images
[est_outliers_uncorr, est_outliers_corr, outlier_tables] = outliers(obj, 'fd', mvmt_matrix);
sum(est_outliers_uncorr) % this shows how many outliers
covs = [covs outlier_tables.outlier_regressor_matrix_uncorr];
% Mean CSF
[gwcsfvalues, gwcsfcomponents] = extract_gray_white_csf(obj);
meancsf = gwcsfvalues(:, 3);
covs = [covs meancsf];
% Here the first component is a head-shift artifact; this will not always
% be the case. You could use PCA or ICA denoising strategies but be
% careful. This code applies to this sample data only.
[component_scores, eigenmap_obj, explained] = pca(obj, 'k', 3);
covs = [covs component_scores(:, 1)];
% High-pass filter - assume TR = 2 sec
[S,KL,KH] = use_spm_filter(2, size(obj.dat, 2),'none', 'specify', 128);
covs = [covs KH];
% Regress all covs out from each voxel and return residuals
% We can ignore colinearity/leverage warnings
obj.X = covs;
regress_output = regress(obj, 'residual', 'noverbose');
obj_denoised = regress_output.resid;
Step 1. Set Up Your Environment
- Download and Install CANLAB Toolboxes:
- Visit CANLAB.github.io and download the latest version of the CANLAB tools.
- Add the CANLAB toolbox directories to your MATLAB (or compatible) path.
- Make sure you have other required dependencies (e.g., SPM if using its filtering routines).
- Prepare Your Inputs:
- Time Series Data: A 4-D NIfTI file per participant (e.g., from fmriprep).
- Anatomical/Segmentation Outputs: CSF, gray matter, and white matter masks (can be extracted using extract_gray_white_csf_tissue_compartments.m).
- Motion and Outlier Information: fmriprep output files or similar that contain motion parameters, displacement measures, and spike regressors.
Step 2. Create the fmri_data Object
The fmri_data object in CANLAB is a container that holds your fMRI data along with metadata and processing options.
Sample data
Orientation to Pinel data files
The Pinel Localizer task was designed to probe several different types of basic cognitive processes, such as visual perception, finger tapping, language, and math. Several of the tasks are cued by reading text on the screen (i.e., visual modality) and also by hearing auditory instructions (i.e., auditory modality). The trials are randomized across conditions and have been optimized to maximize efficiency for a rapid event related design. There are 100 trials in total over a 5-minute scanning session. Read the original paper for more specific details about the task and the dataset paper.
This dataset is well suited for these tutorials as it is (a) publicly available to anyone in the world, (b) relatively small (only about 5min), and (c) provides many options to create different types of contrasts. There are a total of 94 subjects available.
Our Pinel_data_sample_subject is from an individual collected at Dartmouth's Brain Imaging Center (DBIC). It is organized in a particular standardized file format called BIDS (Brain imaging Data structure). The whole dataset is available in BIDS format from Luke Chang's g-node GIN repository.
To have a look at this, let's go to the Pinel_data_sample_subject folder. Drag and drop it into your Command Window, or drag the folder into Terminal (Mac) and copy the path, then use the cd( ) command in Matlab. e.g.,
% This works for me, but you can't copy it exactly unless you have my laptop! Which i hope is not the case.
cd('/Users/torwager/Dropbox (Dartmouth College)/COURSES/Courses_Dartmouth/2021_3_Spring_fMRI_Class/Labs_and_Lectures_for_Instructors/Labs_and_assignments/Lab_5/Pinel_data_sample_subject')
Now, let's look at the files:
ls % list files
The main files we're interested in are:
sub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii
The image data in 4-D Nifti format
sub-sid001567_task-pinel_acq-s1p2_run-03_events.tsv
The task events + their onset times
The rp_.txt file is only created after realignment of the functional images, a preprocessing step. We'll assume this is done already for now (which it was for this sample participant).
After realignment with SPM software (or similar), some new files are created:
rp_sub-sid001567_task-pinel_acq-s1p2_run-03_bold.txt
The "realignment parameters" (rp), head motion estimates, with translation and rotation for each image required to match it to the reference image (often the first image in the series or an average image)
sub-sid001567_task-pinel_acq-s1p2_run-03_bold.mat
Matlab transformation matrices describing the voxel sizes, origin (zero-coordinate), and rotation for each image
rsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii
The realigned image data, with the transformations (i.e., rotation/translation) applied
swrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz
Smoothed (s), warped to MNI template (w), and realigned (r) 4-D image series. This is the one we are going to load and use for this example.
To construct the design matrix for first-level analysis, we're mainly interested in two files: The task (.tsv) file that tells us how to construct the regressors of interest, and the rp*.txt file that we'll use as a set of nuisance covariates.
For further reading and information on file specifications, see the docs on fmriprep preprocessing, or this or this.
% Load the key 4-D image file
fname = which('swrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz');
obj = fmri_data(fname);
Step 3. Build the Nuisance Covariate Model
The goal is to create a design matrix that captures nuisance variance from multiple sources. This matrix typically includes:
1. Intercepts:Use intercept_model.m to generate indicator matrices for each run or session.
This example dataset only includes one run, so we will omit it here, but it is important to include for multi-run data.
2. Motion Parameters: Extract the 24 motion parameters per run (e.g., translations, rotations, and their derivatives/squares) from your fmriprep or SPM outputs.
% 24 movement-related covariates per run
mvmtfname = which('rp_sub-sid001567_task-pinel_acq-s1p2_run-03_bold.txt');
[mvmt_matrix, mvmt_regs_24] = movement_regressors(mvmtfname);
covs = mvmt_regs_24;
Let's plot them:
figure; plot(mvmt_matrix);
legend({'x' 'y' 'z' 'roll' 'pitch' 'yaw'});
title('Movement regressors')
3. Mean CSF Signal: Use either the fmriprep output or the CANLAB function:FD = FD.extract_gray_white_csf();
% Mean CSF
[gwcsfvalues, gwcsfcomponents] = extract_gray_white_csf(obj);
meancsf = gwcsfvalues(:, 3);
covs = [covs meancsf];
This extracts average signals from CSF (and optionally gray and white matter).
4. Spike Regressors:Identify and flag time points with unusual signal spikes (e.g., using FD.outliers).
5. High Displacement Time Points:Incorporate additional covariates for time points with high head displacement, as measured by geometric displacement metrics from fmriprep.
% 'spike' outliers and high framewise-displacement (fd) images
[est_outliers_uncorr, est_outliers_corr, outlier_tables] = outliers(obj, 'fd', mvmt_matrix);
sum(est_outliers_uncorr) % this shows how many outliers
% each outlier or high-fd image gets its own regressor
covs = [covs outlier_tables.outlier_regressor_matrix_uncorr];
6. High-Pass Filter (HPF) Regressors:Generate HPF regressors to remove low-frequency drift.
For example, using SPM’s filtering functions:
cutoff = 128; % Example cutoff in seconds; adjust as needed
[HPreg, ~] = use_spm_filter(FD.dat, FD.TR, cutoff);
Note that we need to know the TR for this! It's 2 sec. In BIDS formatted data, this can be found in the ".json sidecar" file, in the "RepetitionTime" field. We also need the scan length in sec, which in this case is 300 sec or 5 min.
% High-pass filter - assume TR = 2 sec
[S,KL,KH] = use_spm_filter(2, size(obj.dat, 2), 'none', 'specify', 128);
covs = [covs KH];
covs: Combined nuisance regressors into nuisance design matrix X
The various sets of regressors are concatenated into one matrix, covs.
This includes intercept, motion, CSF, spikes, high displacement (movement), HPF
Step 4. Regress Out Nuisance Signals
Now, remove the nuisance variance from each voxel’s time series. The fmri_data.regress method with the 'residual' option performs this operation:
% Regress nuisance signals; output is the residual time series (denoised)
obj_denoised = obj.regress(nuisanceMat, 'residual');
% Regress all covs out from each voxel and return residuals
% We can ignore colinearity/leverage warnings
obj.X = covs;
regress_output = regress(obj, 'residual', 'noverbose');
obj_denoised = regress_output.resid;
After this step, obj_denoised holds the denoised 4-D time series for further analysis.
% Optionally, save the design matrix alongside for documentation
save('nuisance_covariates.mat', 'covs');
% Visualize the cleaned data
plot(obj_denoised)
Step 5. Save the Preprocessed Data
Store the denoised fmri_data object and any relevant nuisance matrices for subsequent analyses:
% Save the fmri_data object
obj_denoised.save('participant_denoised.mat');
This file can be loaded later for connectivity analyses or model fitting.
Step 6. Optional: Explore how nuisance covs overlap with task design
Create a task design matrix from the Pinel task file
Load the realignment parameters and create a set of nuisance covariates
% Load the onsets file into a table object in matlab
taskdata = readtable('sub-sid001567_task-pinel_acq-s1p2_run-03_events.tsv', 'FileType', 'text');
taskdata
% Use a CANlab function to extract the conditions (event types) and which
% events belong to each type
[indic, names] = string2indicator(taskdata.trial_type);
indic = logical(indic);
% Get a cell array of onsets for each condition
[ons, ons_with_durations] = deal(cell(1, length(names)));
for i = 1:length(names)
ons{i} = taskdata.onset(indic(:, i));
ons_with_durations{i} = taskdata.onset(indic(:, i));
ons_with_durations{i}(:, 2) = taskdata.duration(indic(:, i));
end
Now, we've extracted the information in a format we can use to construct a design matrix. Let's try it.
Note that we need to know the TR for this! It's 2 sec. In BIDS formatted data, this can be found in the ".json sidecar" file, in the "RepetitionTime" field. We also need the scan length in sec, which in this case is 300 sec or 5 min.
TR = 2;
X = onsets2fmridesign(ons, TR, 300);
% let's plot it
plotDesign(ons, [], TR);
This returns the task design matrix, X.
We can look at the Variance Inflation Factors to see whether there are mulitcolinearity issues or not. We'll use the 'plot' option and save the handle to this figure so we can re-activate it and add things later:
fig_han = create_figure('vifs'); % Get handle to current figure
vifs = getvif(X, 0, 'plot');
title('Var. inflation factors without nuisance covs')
Now let's see how VIFs are affected by nuisance covs:
fig_han = create_figure('vifs_with_covs'); % Get handle to current figure
vifs = getvif([X covs], 0, 'plot');
% We are interested in the task regs only, so zoom in on those:
set(gca, 'XLim', [0 10], 'YLim', [0 max(10, max(vifs(1:10)) + 0.5)]);
title('Var. inflation factors with nuisance covs')
Step 6. Analysis: Connectivity and Custom RL Models
Once you have the cleaned data, you can proceed with various analyses:
A. Connectivity Analysis
1. Extract Regional Time Series:Use functions such as extract_data (or the CANLAB-supported brainpathway tools) to get time series from regions of interest (ROIs).% Example: extract time series for a specific ROI
roiMask = 'path/to/roi_mask.nii';
roiTS = obj_denoised.extract_data(roiMask);
2. Compute Connectivity Metrics:You can calculate correlation matrices, graph metrics, or other connectivity measures from the extracted time series.
See canlab_connectivity_predict()
B. Custom Reinforcement Learning (RL) Models
1. Prepare the Data for Modeling:Use the denoised voxel or ROI time series as inputs to your RL model.For instance, you might extract trial-by-trial signals aligned to task events.
2. Fit Your RL Model:Implement your RL model (e.g., using MATLAB’s optimization or Bayesian tools). You may code the model in MATLAB or interface with Python (e.g., via MATLAB’s Engine API).
% Pseudocode example for model fitting
modelParams = fitRLModel(roiTS, taskEvents);
3. Evaluate Model Performance:Analyze the model parameters, compare with behavioral data, or run statistical tests (e.g., using fitlme, igls, or CANLAB’s fmri_data.regress for second-level analysis).
Summary
- Data Loading:Read the 4-D NIfTI into an fmri_data object.
- Nuisance Modeling:Build a nuisance regressor matrix using intercepts, motion, CSF, spike, high displacement, and HPF regressors.
- Denoising:Regress out the nuisance signals to create a cleaned dataset.
- Saving:Save the denoised data and nuisance matrices.
- Analyses:
- Connectivity: Extract ROI time series and compute connectivity measures.
- RL Modeling: Fit your custom reinforcement learning model to the cleaned time series.
By following these steps and using the CANLAB toolset, you can efficiently pre-process your fMRI data and set the stage for advanced connectivity and model-based analyses. Happy analyzing!