rs fMRI
After you have carried out your pre-processing and denoising, or, if you are trying to decided which denoising strategy will provide the best quality functional connectivity data, you can run the steps given below. And after have run MRI-QC on your data, and excluded any subjects that have large artefacts, we can dig a little deeper and look at more subtle sources of noise, as well as the impact our chosen preprocessing pipeline, and any additional denoising strategies such as ICA, Tedana or Global Signal Regression (GSR).
If you have run fmriprep on your data, you should have a file that contains a variable known as framewise displacement (FD). For each subject, you should have a FD timeseries equal to the number of volumes. FD is usually defined as the sum of the absolute temporal derivatives of the six motion parameters, following conversion of rotational parameters to distances by computing the arc length displacement on the surface of a sphere with radius 50 mm or above (FD_power). The formula for calculating FD is below (you usually don't need to worry too much about this, as fmriprep has already calculated this for you):
where d denotes translation distances {x,y,z}, and r denotes rotation angles {α,β,γ}. For each participant, a single (scalar) estimate of overall motion, the mean FD, can be calculated by simple taking the average of the FD time series. This variable index the amount of head motion during a scan, and some rules can be used to exclude scans which have above a certain threshold of movement. Note there are slightly different ways to calculate mean FD (FD_Power vs. FD_Jenk) which can result in values that are highly correlated but differ in magnitude. Parkes et al (2019) outline two sets of criteria based on Jenk's method of calculating FD.
- Lenient
- Any scan which has a mean FD (mFD) > 0.55mm is excluded
- Stringent (scans are excluded if any of the below criteria are met)
- mean FD > 0.25 mm
- more than 20% of the FDs were above 0.2 mm
- if any FDs were greater than 5 mm
These are guidelines and can be modified according to your sample. For instance, if you are working in a very high head motion sample, such as in paediatric populations, you may want to be more lenient so you don't have to exclude a large number of scans.
The motion parameters for multiband data is corrupted with respiratory frequencies as highlighted in Power et al. (2019) and Fair et al. (2020).
To filter out these respiratory frequencies a bandpass filter between 0.31hz-0.43hz, this can be implemented with the following script:
function fd = mband_powerFD(mov, TR, head)
% This function computes framewise displacement according to the Power
% method. This code is edited for multiband data to filter out respiratory frequencies according to Fair et al.
% 2020. Movement parameters can be taken from mcflirts .par output
% This script is adapted from GetFDPower.m (Linden Parkes, Alex Fornito.
% Brain & Mental Health Laboratory, 2016) and GetFDJenk_multiband.m (Linden
% Parkes, Tribikram Thapa, Ben Fulcher & Alex Fornito, 2020).
% Updated script:
% Kane Pavlovich, & Alex Fornito, 2022.
% INPUTS
% ------
% mov - an N x 6 matrix containing 6 movement parameters and where
% N = length of the time series.
% TR - TR of scan in seconds
% head - head radius (in mm). Default = 50mm
%
% -------
% OUTPUTS
% -------
% fd - N-1 length vector representing the total framewise
% displacement. There will be k zeros at the start of the
% vector.
if nargin < 3
head=50;
% convert degrees to radians to mm
mov(:,4:6) = head*pi/180*mov(:,4:6);
% % Filter out respiratory frequencies based on Fair et al. 2020
stopband=[0.31 0.41];
nyq = (1/TR)/2; % nyquist
Wn = stopband/nyq; % cutoff frequencies
order = 10; % controls steepness of filter roll-off. Set to same value as in Power et al.2019
[B, A] = butter(order,Wn,'stop'); % construct filter
% detrend and filter motion parameters
mov = detrend(mov,'linear'); % HCP already detrended so no need.
mov = filtfilt(B, A, mov);
% differentiate movement parameters
delta_mov = [
zeros(1,size(mov,2));
diff(mov);
];
% compute total framewise displacement
fd = sum(abs(delta_mov),2);
endNote that Power also suggests downsampling the motion parameters to account for the low TR of multiband data, but this isn't a widely applied step and so we are leaving this out for now until further research on this is conducted.
Carpet plots are becoming an increasingly popular way to visualize subject-level fMRI scan quality. The plot visualizes a matrix of color-coded signal intensities, in which rows represent voxels and columns represent time. The order in which these voxels are displayed has a major effect on the visual interpretation of carpet plots. Please see these two papers for more details:
- Identifying and removing widespread signal deflections from fMRI data: Rethinking the global signal regression problem
- A simple but useful way to assess fMRI scan qualities
What am I looking for?
What you are looking for is large sections of high (white) or low (black) values, known as wide-spread signal deflections (WSD), more details on what causes WSDs can be found here. These are assumed to be indicators of noise and need to be removed or controlled for where possible. Here is an example of a scan where there are large WSDs:
Here is an example where WSDs are not obvious at all:
It is good to look at carpet plots before and after you have applied any preprocessing, to see if your chosen preprocessing strategy had reduced WSDs. An important thing to think about when making and looking at carpet plots is how the voxels are ordered along the Y-axis. Most of the time, these voxels are ordered randomly. You can also order these voxels based on correlation strength to the global signal, or using a clustering technique. Looking at these different orderings can reveal WSDs that are not always visible when viewing carpet plots with random ordering. To see how these ordering methods impact carpet plots, have explore on this page made by Kevin Aquino.
Here is an example of carpet plots before and after some preprocessing and denoising have been applied. Cluster ordering of voxels has been used here to highlight WSDs:
Note that it is still unclear which WSDs represent noise and which do not. Some are clearly tied to head motion whereas others have unknown origins. This paper by the HCP team suggests that some may be of neurophysiological interest: https://pubmed.ncbi.nlm.nih.gov/29753843/ And see also the subsequent debate: https://pubmed.ncbi.nlm.nih.gov/30849529/ https://pubmed.ncbi.nlm.nih.gov/31026516/
So, don't assume that all WSDs should be removed for adequately clean data. If using a technique such as DiCER, one strategy might involve being as unaggressive as possible, while ensuring that QC metrics (discussed below) are on par or better than GSR.
There are some tools already available for making carpet plots:
- CarpetplotR - Bash command line-based tool. Allows for random, global signal and cluster ordering.
Simple download and usage of carpetplotR (see git repo above for more details and trouble shooting):
wget https://raw.githubusercontent.com/sidchop/carpetplotR/main/carpetplotR.R
Rscript carpetplotR.R -f fmri_file.nii.gz
See also:
Sometimes, due to the acquisition parameters of functional MRI scans, certain parts of the the brain do not have adequate signal coverage. Often these are ares of the orbito-frontal cortex, temporal pole or cerebellum. Here is an example of an fMRI image with low signal in frontal and temporal areas:
As a side note, It is not obvious that low-signal == poor signal, but it is seems reasonable to assume that the signal from these areas may contain more noise, and as such, it may be best to exclude these regions. There are a few different ways to quantify which ROIS to exclude.
One way is to calculate each region’s mean BOLD signal across the all scans. Sort these values from largest to smallest, and identify the ‘elbow’ of this distribution using the pairwise differences. This method is outline in Brown et al. 2019 (see supplement). Here is some example R code to implement this method:
Here is an example out put from the code above, where the red line signifies the largest drop in signal, and the regions below this line may be exclude:
Based on the plot above, here is a visualisation of the ROIs that were excluded, overlaid on the mean fmri scan, across all subjects:
Code for a similar method in MatLab can be found here. This is very similar to the Brown et. al method. We also find the 'elbow' of the distribution and use it as the cut-off point, but this is done separatelly for every subject, instead of across all the scans.
In addition to examining scan quality, before and after preprocessing it is also good to visualize the final derivatives on which you will be running your analyses. For instance, if you are doing whole-brain atlas-based connectivity analysis, you should plot the connectivity matrices for each scan, before and after denoising.
Once you have extracted your timeseries, you can read them into your analysis software of choice for further analysis. As we mentioned at the start of this section, functional connectivity matrices (FC matrices; sometimes known as functional connectomes) are often calculated from time series (by calculating the self-correlation) for further analysis. FC matrices can also be good indicators of scan quality.
For example, here is the functional connectivity matrix of a scan before and after GSR has been applied:
What you are looking for here will depend on the type of analysis you are running. For instance, if you are doing whole-brain atlas-based connectivity analysis, you should expect to see a higher correlation within resting-state networks. While there is no ground truth here, we should really see some structure in these matrices, like the ones shown above. If there is no or very little structure, such as the matrix shown below, then there might be an issue with scan quality, or the denoising used.
Likewise, if you are seeing very high correlations across most regions, you might also want to investigate further.
Sometimes, you may find that you can see noise within the correlation structure of these matrices, which might imply that your denoising procedures have not worked correctly, or that you have re-introduced noise into your scan at some point.
The trick is to look at the theory and what is known about the data you are examining. If you are looking at seed-based correlations, do the subject-level seed-to-brain correlation match with what is know about anatomy and function in those networks? Likewise, if you're examining whole-brain resting-state networks connectivity, can you see a resting-state network correlation structure within each subjects data? These are difficult questions, as we do not have ground-truth as to what correctly denoised FC matrices should look like, but I hope the above examples have given you a few guidelines as to what they should not look like.
It can also be useful to look at an average functional connectivity matrix, across all subjects. This is especially useful for comparing the impact of different preprocessing and denoising methods.
A very useful set of metrics for assessing fMRI data quality are known as QC-FC metrics. These metrics look at the correlations between measures of quality (such as mean Frame-wise Displacement) and functional connectivity. Ideally, there should be a low correlation between mean Frame-wise Displacement (a proxy for head motion) and functional connectivity values. For a more detailed discussion on QC-FC metrics, see Parkes et al (2019) and Cric et al (2018) - the Ciric paper come with some code you can use as well, although, I provide similar R code below.
FD-FC Correlation First, we can assess the correlation between each subjects connectivity value (Pearson's r) and mean FD, at each edge. Here is some R code which reads in a vector of mean FD values, and each subjects connectivity matrix, and calculates this FD-FC correlation, and outputs a .png file.
FD-FC-Distance Correlation We can also examine the extent to which FD–FC correlations are dependent on the Euclidean distance between the centroids of any pair of regions. The presence of head-motion artificially increases correlations between regions that are closer, and decreases correlation between regions that are further part. Therefore, if there is no or minimal impact of head-motion on FC, we should not see any pattern of correlation between distance and FD-FC correlations. For instance, have a look at this FD-FC-Distance plot from Parkers et al:
As we can see, there is a clear trend for higher correlations between FD and FC at shorter distances, and lower correlations at longer distances. Now here is an example where a particular denoising strategy has substantially reduced this association (which is what we hope to see):
If you do now have the centroid coordinates for the atlas you are using, you can use FSL (fslstats) to get centroid/COG coordinates from .nii atlas or use the R code below:
library(neurobase)
atlas <- readNIfTI2('path/to/atlas.nii')
message('Calulating center of gravity (COG) for each ROI in the atlas image.')
#get voxel coordinates
cog <- matrix(nrow=length(unique(atlas[atlas>0])), ncol = 3)
for (s in 1:length(unique(atlas[atlas>0]))){
temp <- atlas
temp[temp!=s] <- NA
cog[s,] <- neurobase::cog(temp,)
}
#vox to MNI
for (c in 1:dim(cog)[1]){
cog[c,] <- oro.nifti::translateCoordinate(cog[c,],nim=atlas)
}
cog <- round(cog, 0)
write.table(cog, paste0('/atlas_centroids.txt'),
col.names = F, quote = F, row.names = F)
Below is some R code to make FD-FC-Distance plots similar to the ones above:
library(ggplot2)
library(bio3d)
dist <- read.table(opt$cog, header = F) #load centroid coords
xyz.dist <- dist.xyz(dist) #calculate euclidian distance between each pair of centroids
distance <- xyz.dist[upper.tri(xyz.dist)]
fc_fd <- read.table(opt$fdfc, header = F)
fc_fd <- fc_fd[upper.tri(fc_fd)]
xy <- cbind(distance,fc_fd)
message('Making FC-FC distance plots')
ggplot(as.data.frame(xy), aes(x=distance, y=fc_fd)) +
geom_point(alpha=0.3)+
# ggtitle("fix_gsr") +
stat_summary_bin(fun="mean", bins=20,
color='orange', size=2, geom='point') +
stat_summary_bin(fun.data = mean_sdl, bins=20,
color='orange', size=0.5) + ylim(-1,1)
ggsave(filename = paste0(opt$out,'/group/FD_FC_dist.png'),
device = "png",
width = 7,
height = 5)
Example output:
A useful final check to conduct is to examine the subject to subject correlation, between the vectorised FC values for each subject. This allows us to identify subjects/scans where FC values do not correlate with other subjects/scans, or potentially highly correlated clusters of subjects. Here is an example of such a subject by subject correlation matrix for s study with 80 subjects:
We can see that there is at-least one subject where their FC values are weakly correlated with all other subjects. It might be worth examining this subject more closely (e.g. look at the FC matrix, carpet plot, FD values). It can also be informative to make these plots before and after applying your denoising strategies. For instance, here is the same sample from above, this time after some more denoising:
As we can see, the subject no longer 'sticks out' and most subjects are moderately correlated.