% GLM Analysis Script for fMRI Data using SPM
%
% This script performs a General Linear Model (GLM) analysis on fMRI data that has been preprocessed using fMRIprep and
% is organized in BIDS format. It is designed to run in MATLAB using SPM, specifically for researchers who may be
% unfamiliar with SPM scripting but have fMRI data analysis needs. The script follows standard steps for GLM setup,
% specification, estimation, and contrast definition, which are essential in understanding task-related brain activity.
%
% Overview of What This Script Does:
% - Sets up paths and parameters.
% - Loads event and confound data for each run and prepares it for GLM analysis.
% - Specifies the GLM parameters, including timing, high-pass filter, and basis function.
% - Runs model specification and estimation to fit the GLM to the fMRI data.
% - Defines contrasts to test specific hypotheses.
% - Saves the final model and a copy of this script for future reference.
%
% Steps in this Script:
% 1. **Set Parameters**: Define paths to data, fMRI timing parameters, selected tasks, contrasts, and subjects.
% 2. **Load Events and Confounds**: Each run`s event (trial timing) and confound (noise regressors) data are loaded
% from BIDS-compliant TSV and JSON files produced by fMRIprep. This prepares each run for GLM analysis.
% 3. **SPM Model Specification**: Sets up a structure (called`matlabbatch`) that stores all necessary parameters
% for defining the GLM model. Each parameter here mirrors options in the SPM GUI.
% 4. **Run Model Specification and Estimation**: Once all parameters are set in`matlabbatch`, we run batches 1
% and 2. Batch 1 specifies the model, and batch 2 estimates the model using the data and settings we’ve defined.
% 5. **Check SPM.mat File**: Verifies the existence of the SPM.mat file after running the model, confirming the
% model specification and estimation steps were successful.
% 6. **Define and Run Contrasts**: Sets up contrasts based on hypotheses about brain activity, applies them, and
% stores results. This is done in batch 3.
% 7. **Save Script**: Saves a copy of this script in the output folder for reproducibility and future reference.
% Working with matlabbatch and SPM Batches:
% - `matlabbatch` is a structure used by SPM to store each step of an analysis, organized as separate fields for
% each processing stage (e.g., model specification, model estimation, contrast definition). Each entry in
% `matlabbatch` corresponds to a batch operation you could run manually in the SPM GUI. SPM’s batch system is
% powerful because it allows us to set parameters programmatically in a way that is identical to the GUI.
% - **Setting Parameters**: Each field in `matlabbatch` represents a setting or action, matching options in the
% SPM GUI. For example, in this script, `matlabbatch{1}.spm.stats.fmri_spec.timing.RT` corresponds to setting
% the Repetition Time (RT) in the SPM GUI.
% - **Retrieving GUI Code**: To see the code for any parameter set in the SPM GUI, configure the GUI settings and
% then select **File > Save Batch and Script**. SPM will generate MATLAB code that mirrors the GUI settings,
% which you can use to understand how `matlabbatch` fields align with GUI options.
%
% Understanding SPM Batches and Sessions:
% - In SPM, **sessions** represent individual runs in a study. Each run (session) is specified separately within
% `matlabbatch`, which enables you to specify different conditions or timing settings per run if necessary.
% - We use indexing within `matlabbatch` to indicate which session or contrast we’re defining. For instance,
% `matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).scans` refers to the scans (images) in session `runIdx`.
%
% How SPM Executes Batches:
% - SPM uses a deferred execution approach: parameters are specified first by configuring `matlabbatch`, and
% then `spm_jobman('run', matlabbatch)` initiates all specified processing steps in sequence, similar to
% clicking "Run" in the GUI. In this script, we define `matlabbatch` elements in sequence (batches 1 to 3)
% and run them together in one command to automate the workflow.
%
% Assumptions about Input Data:
% - Data must be preprocessed using fMRIprep and follow the BIDS format.
% - Event files (in TSV format) must adhere to BIDS specifications.
% - Confound files (in JSON and TSV formats) must adhere to fMRIprep output standards.
%
% Requirements:
% - SPM12 (or a compatible version) must be installed and added to the MATLAB path.
%
% This script automates the process typically done in the SPM GUI.
% Parameters:
% - selectedSubjectsList: List of subject IDs to include in the analysis. Can be a list of integers or '*' to include all subjects.
% - selectedRuns: List of run numbers to include in the analysis. Can be a list of integers or '*' to include all runs.
% - selectedTasks: Structure array with information about the task(s) to analyze.
% Each task must have the following fields:
% - name: String. The name of the task.
% - contrasts: Cell array of strings. The name(s) of the contrast(s).
% - weights: Cell array of structs. The weights for the contrast(s).
% - smoothBool: Boolean. Whether to smooth the data before the GLM.
%
% Example of defining tasks, weights, and contrasts:
% selectedTasks(1).name = 'loc';
% selectedTasks(1).contrasts = {'Faces > Objects', 'Objects > Scrambled', 'Scenes > Objects'};
% selectedTasks(1).weights(1) = struct('Faces', 1, 'Objects', -1, 'Scrambled', 0, 'Scenes', 0);
% selectedTasks(1).weights(2) = struct('Faces', 0, 'Objects', 1, 'Scrambled', -1, 'Scenes', 0);
% selectedTasks(1).weights(3) = struct('Faces', 0, 'Objects', -1, 'Scrambled', 0, 'Scenes', 1);
% selectedTasks(1).smoothBool = true;
%
% Paths:
% - fmriprepRoot: Path to the root folder of the fMRIprep output.
% - BIDSRoot: Path to the root folder of the BIDS dataset.
% - outRoot: Path to the output root folder for the analysis results.
%
% Preprocessing Options (uses fMRIprep confounds table):
% - pipeline: Denoising pipeline strategy for SPM confound regression.
% It should be a cell array of strings specifying the strategies to use.
% Possible strategies:
% HMP - Head motion parameters (6,12,24)
% GS - Brain mask global signal (1,2,4)
% CSF_WM - CSF and WM masks global signal (2,4,8)
% aCompCor - aCompCor (10,50)
% MotionOutlier - motion outliers FD > 0.5, DVARS > 1.5, non-steady volumes
% Cosine - Discrete cosine-basis regressors (low frequencies) -> HPF
% FD - Raw framewise displacement
% Null - Returns a blank data frame
%
% Example:
% pipeline = {'HMP-6','GS-1','FD'};
% This selects the 6 head motion parameters, the global signal, and the raw Framewise Displacement value.
%
% Author: Andrea Ivan Costantino
% Date: 5 July 2023
% This section of the script defines parameters and settings necessary for a GLM analysis in SPM, using fMRI data
% preprocessed with fMRIPrep. We first set up the required paths, select specific subjects, tasks, and contrasts,
% and define analysis parameters like repetition time and timing information. The final block iterates through
% subjects and tasks, setting up each run of the GLM and checking for the existence of the required data before
% beginning the analysis. This setup ensures a consistent pipeline for different subjects and tasks.
clc
clear
%% PARAMETERS
% Define the space of the images ('T1w' or 'MNI'); in this case, we use MNI space at 2mm resolution.
niftiSpace = 'MNI152NLin2009cAsym_res-2';
% Define the root directory for the BIDS dataset (where the data is stored).
BIDSRoot = '/data/projects/chess/data/BIDS';
% Define paths relative to the BIDS root:
fmriprepRoot = fullfile(BIDSRoot, 'derivatives', 'fmriprep'); % Path to fMRIprep output data.
outRoot = fullfile(BIDSRoot, 'derivatives', 'fmriprep-SPM-TEST', niftiSpace); % Directory for storing analysis results.
tempDir = fullfile(BIDSRoot, 'derivatives', 'fmriprep-preSPM'); % Temporary storage for intermediate files.
% Define fMRI timing parameters.
RT = 2; % Repetition time (TR) in seconds, the time between consecutive volume acquisitions.
fmri_t = 60; % Microtime resolution (number of slices).
fmri_t0 = fmri_t / 2; % Microtime onset (reference slice for alignment).
% Define the subjects and runs to include in the analysis.
selectedSubjectsList = [1,2,3]; % List of subject IDs to include (or use '*' for all subjects).
selectedRuns = '*'; % Runs to include, with '*' indicating all available runs.
% Define tasks with contrasts and associated weights for GLM analysis.
% Each task contains fields: name, contrasts, weights, and smoothBool.
% Task 1: Localizer Task
selectedTasks(1).name = 'loc'; % Name of task as used in BIDS.
selectedTasks(1).contrasts = {'Faces > Objects', 'Objects > Scrambled', 'Scenes > Objects'}; % Contrasts of interest.
selectedTasks(1).weights(1) = struct('Faces', 1, 'Objects', -1, 'Scrambled', 0, 'Scenes', 0); % Weights for contrast 1.
selectedTasks(1).weights(2) = struct('Faces', 0, 'Objects', 1, 'Scrambled', -1, 'Scenes', 0); % Weights for contrast 2.
selectedTasks(1).weights(3) = struct('Faces', 0, 'Objects', -1, 'Scrambled', 0, 'Scenes', 1); % Weights for contrast 3.
selectedTasks(1).smoothBool = true; % Apply smoothing before GLM for this task.
% Task 2: Experimental Task
selectedTasks(2).name = 'exp';
selectedTasks(2).contrasts = {'Check > No-Check'};
selectedTasks(2).weights(1) = struct('check', 1, 'nocheck', -1); % Weights for contrast 1.
selectedTasks(2).smoothBool = false; % Do not apply smoothing for this task.
% Specify confound regressors to include in the GLM.
pipeline = {'HMP-6','GS-1','FD'}; % Use head motion parameters, global signal, and framewise displacement.
% Find folders for each subject based on the selected subjects list.
sub_paths = findSubjectsFolders(fmriprepRoot, selectedSubjectsList);
%% RUN THE GLM FOR EACH SUBJECT AND TASK
% Loop through each subject in the list.
for sub_path_idx = 1:length(sub_paths)
% Loop through each defined task for analysis.
for selected_task_idx = 1:length(selectedTasks)
% Clear previous configurations in matlabbatch to avoid conflicts.
clearvars matlabbatch
%% SUBJECT AND TASK INFO
% Get the current task details (name, contrasts, smoothing preference).
selectedTask = selectedTasks(selected_task_idx).name; % Task name.
contrasts = selectedTasks(selected_task_idx).contrasts; % Task contrasts.
smoothBool = selectedTasks(selected_task_idx).smoothBool; % Boolean: smooth images before GLM?
% Extract the subject ID from the folder name (formatted as 'sub-01').
subPath = sub_paths(sub_path_idx);
subName = subPath.name; % e.g., 'sub-01'
sub_id = strsplit(subName,'-'); % Splits 'sub-01' into {'sub', '01'}
selectedSub = str2double(sub_id{2}); % Extracts and converts the subject ID to a number.
% Define the output path for storing GLM results for this subject and task.
outPath = fullfile(outRoot, 'GLM', subName, selectedTask);
% Define paths to the functional data for this subject.
funcPathSub = fullfile(fmriprepRoot, subName, 'func'); % Functional data path.
bidsPathSub = fullfile(BIDSRoot, subName, 'func'); % BIDS formatted path.
% Check if the specified task exists in the subject’s data.
files = dir(funcPathSub); % List files in the functional data directory.
fileNames = {files.name}; % Extract file names.
containsTask = any(contains(fileNames, ['task-', selectedTask])); % Check for the task in file names.
% If the task is missing for this subject, skip to the next iteration.
if ~containsTask
warning(['Task ', selectedTask, ' not found for ' subName ' in ' funcPathSub '. Skipping...']);
continue;
end
% Print a status update showing the subject and task being processed.
fprintf('############################### \n# STEP: running %s - %s #\n############################### \n', subName, selectedTask);
% This first section finds and loads the event and confound files for each run, based on selections specified earlier.
% We retrieve and organize the event data (which includes onset times and durations) and confound data
% (like motion parameters or physiological noise) for each run of a task. These inputs are essential for setting up the
% GLM model, as they provide information needed for modeling the BOLD signal in relation to task events and controlling
% for noise.
%% FIND AND LOAD EVENTS AND CONFOUNDS FROM FMRIPREP FOLDER
% Retrieve event and confound files for selected runs.
if ismember('*', selectedRuns)
% If '*' is specified, include all available runs for this subject and task.
eventsTsvFiles = dir(fullfile(bidsPathSub, strcat(subName, '_task-', selectedTask, '_run-*_events.tsv')));
json_confounds_files = dir(fullfile(funcPathSub, strcat(subName, '_task-', selectedTask, '_run-*_desc-confounds_timeseries.json')));
tsv_confounds_files = dir(fullfile(funcPathSub, strcat(subName, '_task-', selectedTask, '_run-*_desc-confounds_timeseries.tsv')));
else
% Include only the specific runs listed in `selectedRuns`.
eventsTsvFiles = arrayfun(@(x) dir(fullfile(bidsPathSub, strcat(subName, '_task-', selectedTask, '_run-', sprintf('%01d', x), '_events.tsv'))), selectedRuns, 'UniformOutput', true);
json_confounds_files = arrayfun(@(x) dir(fullfile(funcPathSub, strcat(subName, '_task-', selectedTask, '_run-', sprintf('%01d', x), '_events.json'))), selectedRuns, 'UniformOutput', true);
tsv_confounds_files = arrayfun(@(x) dir(fullfile(funcPathSub, strcat(subName, '_task-', selectedTask, '_run-', sprintf('%01d', x), '_events.tsv'))), selectedRuns, 'UniformOutput', true);
end
% Sort the retrieved files alphabetically by name for consistent ordering.
eventsTsvFiles = table2struct(sortrows(struct2table(eventsTsvFiles), 'name'));
json_confounds_files = table2struct(sortrows(struct2table(json_confounds_files), 'name'));
tsv_confounds_files = table2struct(sortrows(struct2table(tsv_confounds_files), 'name'));
% Ensure the number of event and confound files matches across all run files.
assert(numel(eventsTsvFiles) == numel(json_confounds_files) && numel(json_confounds_files) == numel(tsv_confounds_files), ...
'Mismatch in number of TSV events, TSV confounds, and JSON confounds files in %s', funcPathSub)
% This second section sets non-run-specific SPM model parameters, defining the core structure of the GLM model.
% Here we specify settings such as the repetition time, microtime resolution, basis functions, and whether global
% normalization should be applied. These parameters remain consistent across different runs and subjects, creating
% a standardized foundation for the GLM analysis.
%% SPM MODEL PARAMETERS (NON RUN DEPENDENT)
% Define non-run-specific parameters for the GLM in SPM.
matlabbatch{1}.spm.stats.fmri_spec.dir = cellstr(outPath); % Directory where SPM output will be saved.
matlabbatch{1}.spm.stats.fmri_spec.timing.units = 'secs'; % Units for event timing (set to seconds).
matlabbatch{1}.spm.stats.fmri_spec.timing.RT = RT; % Repetition time (TR) defines time between successive scans.
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t = fmri_t; % Microtime resolution - number of time bins within one TR.
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = fmri_t0; % Microtime onset - reference slice for alignment.
matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {}); % No factorial design applied here.
matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [0 0]; % Canonical HRF without derivatives (time/dispersion).
matlabbatch{1}.spm.stats.fmri_spec.volt = 1; % Volterra interactions disabled (linear model).
matlabbatch{1}.spm.stats.fmri_spec.global = 'None'; % No global normalization applied.
matlabbatch{1}.spm.stats.fmri_spec.mthresh = 0.8; % Masking threshold excludes voxels below 80% of mean signal.
matlabbatch{1}.spm.stats.fmri_spec.mask = {''}; % No explicit mask specified.
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)'; % AR(1) autocorrelation correction for temporal dependencies.
% Set up parameters for estimating the GLM in the second stage.
matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model specification: SPM.mat File', ...
substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
matlabbatch{2}.spm.stats.fmri_est.write_residuals = 0; % Do not save residuals.
matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1; % Use classical estimation method (ReML).
% This section iterates over each run of the task for the current subject, setting up the required SPM parameters
% for each run. We load event and confound data, identify the correct NIfTI files, apply smoothing if required,
% and configure session-specific settings in the SPM model specification. Each step is verified to ensure data
% integrity, and any issues with missing or multiple files are flagged for user review.
%% SPM RUN SETTINGS (E.G., EVENTS, CONFOUNDS, IMAGES)
% Loop over each run of the task.
for runIdx = 1:numel(eventsTsvFiles)
% Load event information into a structure compatible with SPM.
events_struct = eventsBIDS2SPM(fullfile(eventsTsvFiles(runIdx).folder, eventsTsvFiles(runIdx).name));
selectedRun = findRunSubstring(eventsTsvFiles(runIdx).name); % Identify the current run (e.g., 'run-01').
fprintf('## STEP: TASK %s - %s\n', selectedTask, selectedRun)
% Select the corresponding confound files for the current run.
jsonRows = filterRowsBySubstring(json_confounds_files, selectedRun);
confoundsRows = filterRowsBySubstring(tsv_confounds_files, selectedRun);
% Check for a single JSON confounds file, and flag errors if multiple or none are found.
if length(jsonRows) ~= 1
error('More than one JSON file found for the specified run. Please check the dataset.');
elseif isempty(jsonRows)
error('No JSON file found for the specified run. Please check the dataset.');
else
jsonRow = jsonRows{1, :};
end
% Check for a single TSV confounds file, and handle errors similarly.
if length(confoundsRows) ~= 1
error('More than one TSV confounds file found for the specified run. Please check the dataset.');
elseif isempty(confoundsRows)
error('No TSV confounds file found for the specified run. Please check the dataset.');
else
confoundsRow = confoundsRows{1, :};
end
% Extract confound data based on the specified pipeline, loading it into an array.
confounds_array = fMRIprepConfounds2SPM(fullfile(jsonRow.folder, jsonRow.name),...
fullfile(confoundsRow.folder, confoundsRow.name), pipeline);
% Define the pattern to locate the NIfTI file and check for its existence.
filePattern = strcat(subName, '_task-', selectedTask, '_', selectedRun, '_space-', niftiSpace, '_desc-preproc_bold');
niiFileStruct = dir(fullfile(fmriprepRoot, subName, 'func', strcat(filePattern, '.nii')));
% Handle cases with missing or multiple NIfTI files.
if size(niiFileStruct, 1) > 1
error('Multiple NIFTI files found for %s.', selectedRun)
elseif isempty(niiFileStruct)
% If no NIfTI files are found, check for compressed (.nii.gz) files and decompress if necessary.
niiGzFilePattern = fullfile(funcPathSub, strcat(filePattern, '.nii.gz'));
niiGzFileStruct = dir(niiGzFilePattern);
if isempty(niiGzFileStruct)
warning('No NIFTI file found for run %s. SKIPPING!', selectedRun)
continue
elseif size(niiGzFileStruct, 1) > 1
error('Multiple NIFTI.GZ files found for this run.')
else
% Decompress the .nii.gz file if found.
niiGzFileString = fullfile(niiGzFileStruct.folder, niiGzFileStruct.name);
gunzippedNii = gunzipNiftiFile(niiGzFileString, fullfile(tempDir, 'gunzipped', subName));
niiFileStruct = dir(gunzippedNii{1});
end
end
% Get the full path of the (possibly decompressed) NIfTI file.
niiFileString = fullfile(niiFileStruct.folder, niiFileStruct.name);
% Smooth the NIfTI file if smoothing is enabled for this task.
if smoothBool
niiFileString = smoothNiftiFile(niiFileString, fullfile(tempDir, 'smoothed', subName));
else
fprintf('SMOOTH: smoothBool set to false for this task. Skipping...\n')
end
% Specify functional images in SPM model specification for the current run.
niiFileCell = {niiFileString};
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).scans = spm_select('expand', niiFileCell); % Expanded list of images.
% Set high-pass filter (HPF) to a value higher than the run duration to avoid filtering.
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).hpf = (matlabbatch{1}.spm.stats.fmri_spec.timing.RT * ...
size(matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).scans, 1)) + 100;
% Add conditions (events) to the model for the current run.
for cond_id = 1:length(events_struct.names)
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).cond(cond_id).name = events_struct.names{cond_id};
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).cond(cond_id).onset = events_struct.onsets{cond_id};
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).cond(cond_id).duration = events_struct.durations{cond_id};
end
% Add confound regressors to the model for the current run.
for reg_id = 1:size(confounds_array, 2)
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).regress(reg_id).name = confounds_array.Properties.VariableNames{reg_id};
matlabbatch{1}.spm.stats.fmri_spec.sess(runIdx).regress(reg_id).val = confounds_array{:, reg_id};
end
end
% This final section initializes the SPM defaults, executes the GLM model specification and estimation,
% verifies that the model ran successfully by checking for the SPM.mat file, and then defines and
% applies contrasts. Finally, the script is copied to the output folder for replicability. Each step
% is verified with status messages, providing clear feedback for the user.
%% RUN BATCHES 1 AND 2
% Initialize SPM defaults for fMRI analysis.
spm('defaults','fmri'); % Load SPM defaults for fMRI.
spm_jobman('initcfg'); % Initialize the SPM job configuration.
% Execute the model specification (batch 1) and model estimation (batch 2).
fprintf('GLM: Running GLM for: %s - TASK: %s\n', subName, selectedTask)
spm_jobman('run', matlabbatch(1:2));
fprintf('GLM: DONE!\n')
%% FIND SPM.mat
% After running the GLM, check for the existence of the SPM.mat file in the output path.
% The SPM.mat file is essential for contrast specification.
spmMatPath = fullfile(outPath, 'SPM.mat');
if ~exist(spmMatPath, 'file')
error('SPM.mat file not found in the specified output directory.');
end
%% Save Boxcar plot and design matrix of estimated model
SPMstruct = load(spmMatPath);
plotBoxcarAndHRFResponses(SPMstruct, outPath)
saveSPMDesignMatrix(SPMstruct, outPath)
%% RUN BATCH 3 (CONTRASTS)
% Define the path to the SPM.mat file in the batch structure for contrast estimation.
matlabbatch{3}.spm.stats.con.spmmat(1) = {spmMatPath};
% Loop through each contrast defined for the current task.
for k = 1:length(contrasts)
% Adjust contrast weights to fit the current design matrix structure.
weights = adjust_contrasts(spmMatPath, selectedTasks(selected_task_idx).weights(k));
% Define each contrast and its properties in the SPM batch structure.
matlabbatch{3}.spm.stats.con.consess{k}.tcon.weights = weights; % Contrast weights.
matlabbatch{3}.spm.stats.con.consess{k}.tcon.name = contrasts{k}; % Name of the contrast.
matlabbatch{3}.spm.stats.con.consess{k}.tcon.sessrep = 'none'; % No session-specific replication.
end
% Execute the contrast batch (batch 3).
fprintf('Setting contrasts...\n');
spm_jobman('run', matlabbatch(3));
fprintf('DONE!\n');
%% Save Contrast Plots
% Thresholds for statistical analysis
thresholds = {
0.001, ...
0.01, ...
0.05 ...
};
% Iterate over contrasts to generate plots
for constrastIdx = 1:length(selectedTasks(taskIndex).contrasts)
% Iterate over thresholds to generate plots
for thresholdIndex = 1:length(thresholds)
% Generate and Save Contrast Overlay Images
% Set crosshair coordinates (modify if needed)
crossCoords = [40, -52, -18]; % FFA. Change as needed.
% Set the index of the contrast to display (modify if needed)
spmContrastIndex = constrastIdx;
% Call the function to generate and save contrast overlay images
generateContrastOverlayImages(spmMatPath, outPath, fmriprepRoot, subjectName, pipelineStr, thresholds{thresholdIndex}, spmContrastIndex, crossCoords);
end
end
end
%% SAVE SCRIPT
% Save a copy of this script in the output directory for documentation and replicability.
FileNameAndLocation = [mfilename('fullpath')]; % Get full path to this script.
script_outdir = fullfile(outPath,'spmGLMautoContrast.m'); % Set destination path for script copy.
currentfile = strcat(FileNameAndLocation, '.m'); % Get filename with extension.
copyfile(currentfile, script_outdir); % Copy script to output directory.
end