Noise modelling with PhysIO toolbox - GLM masking threshold

[SaraCal]

Dear Lars @mrikasper ,

I contacted you a long time ago to ask for help regarding noise modelling using the PhysIO toolbox. I now have a quick follow-up question.

Back then I was analysing resting-state data, and the way I’ve set the preprocessing was by including the PhysIO toolbox in the pipeline after Normalisation, then having a GLM to regress out the noise, and smoothing the output of that GLM. You suggested to use a masking threshold of 0.5 in this first GLM, in order to avoid getting holes in points of low intensity. I then used a threshold of -Inf in the second GLM (the actual first level GLM after @preprocessing).

I am about to perform a similar preprocessing pipeline on some new data, which are task data this time. I was wondering if you still recommend setting up a conservative threshold of 0.5 in the GLM? Just trying to understand if there are new recommendations (we spoke more than 2 years ago) and if something needs to be different for task data compared to resting-state.

Thank you for your help!

Best wishes, Sara

[mrikasper]

Dear Sara,

very nice to hear from you, and to see that you are still working with PhysIO! I hope everything is well despite the toll that the pandemic has taken on all of us in the last 2 years.

• The point of a liberal masking threshold is still valid, regardless of task- vs resting-state analysis. It becomes more important with higher field strength (e.g., at 7T I recommend thresholds of 0.2 or lower) or with higher-channel head coils (e.g., 64 vs 32 or 16 receive channels), because the intensity variation is larger and can vary a lot around the mean (which is what the threshold refers to, i.e., 50% of the mean pixel intensity in the image).

• What does change however is the order of preprocessing and noise modeling. In a task-based analysis, I would always include the physiological noise regressors as “multiple regressors” into your task-GLM, so have only one 1st level GLM step. This is, because the (anti-)correlations between task and physiological regressors may create false negatives or positives, if you regress them out separately.

  • I have given a talk on physiological noise modeling a while back at ISMRM, where I explain this a bit more: https://www.youtube.com/watch?v=ffq2HS8qs1M

• And in general, if you haven’t done so, please use the most current version of PhysIO (v8.0.1), because we have made quite some improvements, in particular to preprocessing of the physiological recordings and modeling of respiratory volume per time in the last 2 years (see, e.g., this paper).

I hope this helps, good luck with your study!

All the best, Lars

[SaraCal]

Dear Lars, thank you very much for the useful info! I have a quick follow-up question: in my previous (resting-state) pipeline, I was writing the residuals after the first GLM that was included in the preprocessing, and was using them as input to the next GLM (which was still a within-subject GLM).

I will not have to do the same now (i.e. using the residuals as input for the next GLM) since I only have one 1st-level GLM and then a 2nd-level (group) GLM, right?

Thank you for your help!!

Sara

[mrikasper]

Dear Sara,

you are right, everything (task+PhysIO modeling) goes into the same first level GLM.

All the best, Lars

[SaraCal]

Dear Lars, thank you for the confirmation. I tested my pipeline on one pilot participant and I now want to check for the outcomes of denoising using tapas_physio_report_contrasts() and tapas_physio_compute_tsnr_gains().

I have 4 runs of task per participant so I perform the preprocessing on each run separately, and then I compute some contrasts in the 1st-lev GLM while averaging (replicate the contrasts) across the 4 runs. I therefore have one SPM.mat file but 4 physio.mat files from the preprocessing. For this reason I don't understand if I need to use the tapas_physio functions 4 times (for each run) or only once. What do you think?

For now I am performing the check for each run (see my code attached), but I get the following error message when computing the tsnr gains:

Operands to the || and && operators must be convertible to logical scalar
values.
 
Error in tapas_physio_compute_tsnr_spm (line 129)
if isInvertableContrast && doInvert
 
Error in tapas_physio_compute_tsnr_spm (line 188)
        tSnrCompareImage = tapas_physio_compute_tsnr_spm(...
 
Error in tapas_physio_compute_tsnr_spm (line 188)
        tSnrCompareImage = tapas_physio_compute_tsnr_spm(...
 
Error in tapas_physio_compute_tsnr_gains (line 115)
        [tSnrImageArray{c}, fileTsnrArray{c}, ...
 
Error in check_denoising_quality (line 60)
        [tsnrgainarray, tsnrarray] =
        tapas_physio_compute_tsnr_gains([thissessionpath,'/physio.mat'],'SPM.mat');

Finally, the contrast report works fine but I'd just like to ask you whether the quality of the denoising is ok, since the F maps basically show the whole brain...I attached one of the contrast reports too in the zip file but if you need more info I have put my SPM.mat and physio.mat files in a drive folder (along with all the contrast images, anatomical file, etc..):

https://drive.google.com/drive/folders/1DKobdFaTThG4ZaGvUMJCp2H_aOxAiKe2?usp=sharing

Thank you very much for your help!

Sara

code_and_contrast_report.zip

[mrikasper]

Dear Sara,

I will have a look at the data you sent, thank you for preparing it!

1. From a first glance, the F-contrasts for the noise rois indeed look too widespread.
2. For motion regressors, there is often a whole brain effect, but the "corona" at the brain edges looks a bit unfamiliar.
   • I must say, however, I am not so used to multi-run data, maybe it's all explained by the higher degrees of freedom you have.
3. Overall the images also look a bit smooth, can you remind me of the resolution of the raw data and FWHM of any smoothing you used?
4. For the tSNR computation, could you try to replace the && in line 129 of tapas_physio_compute_tsnr_spm by &, please, and see what happens?

All the best, Lars

[SaraCal]

Dear Lars,

thank you for your quick reply. 3. The original resolution is 2.5x2.5x2.5 mm and I used a FWHM of 5x5x5 mm for the smoothing 4. seems like it solved the issue! Currently running and I can see it's generating the tsnr ratio images

P.S. in case it helps: I also tried to denoise the data with slightly different parameters (95% threshold instead of 99% and 3 principal components instead of 5) but I got very very similar results from the contrast report (if not the same)

Thank you for looking into this!

Best, Sara

[SaraCal]

Hi again Lars,

sorry for sending many messages, I just wanted to quickly update you on what I've attempted in case it helps:

I've re-done the GLM but this time for each session separately instead of averaging across sessions, purely just to check whether the multi-run GLM could somehow affect the quality check. Seems to me like it didn't change the outcomes that much - I am attaching the contrast reports.

Best wishes, Sara

physio_contrasts.zip

[SaraCal]

Dear Lars, apologies for writing again - any idea of what could cause the denoising to be weird? Let me know if you need more info from me, thank you so much!

Best wishes, Sara

[mrikasper]

Dear Sara,

sorry for the delay, end-of-year wrap-ups everywhere.

Thank you for sending the individual session contrast reports, I find them much easier to interpret. So the combination of the 4 sessions in the first report might obscur individual differences (e.g., due to different motion magnitude in different sessions), or the high degrees of freedom enabling every regressor to explain significant variance anywhere.

Most of them look reasonable, in particular in session 2 and 3 you see very clear brain edge effects of the motion regressors. Also, for the physiological regressors circle of Willis, ventricles, ACC and straight sinus seem to be high variance regions, as expected.

The single session reports also have the advantage that I see the design matrix properly, so I had two extra questions about it:

1. The first 4 regressors are the task regressors, correct? Looks like an event related design with variable inter-trial-intervals? And you are only modeling the canonical HRF, no derivatives, correct?
2. In the first run, at about volume 160, it looks as if there is a sudden spike in the regressors. Did you have a look at the motion regressors (also in general, for the other runs) and try out, e.g., a higher verbosity level of PhysIO and enabling the "censoring" option? This would give you some insight on the motion magnitude of this subject, I am wondering whether they moved a lot and that might explain some wide-spread variance.

All the best, Lars

[SaraCal]

Dear Lars, many thanks for your reply! Glad to hear the single session reports look better. Regarding your extra questions: 1) you're right that I have an event-related design with variable ITI and I am only modeling the canonical HRF. 2) I am enabling censoring in the physio toolbox using the MAXVAL approach with threshold 2 mm / 2 degrees but I have left the verbosity level as default. I can try increasing that but I have a follow-up question: how can I check for the outcomes of denoising if the multi-run GLM is an issue? With this pilot I re-did the GLM for each run but I wouldn't do that for all the actual participants if there was perhaps another way of checking

Thank you very much!

Best wishes, Sara

[SaraCal]

Dear Lars,

hope you are well! I am performing noise modelling on the same task data (I've collected some more data now but all the parameters are the same as stated in the previous posts of this thread) and I noticed I get the following warning message when the PhysIO toolbox creates the noise regressors:

05-Apr-2022 18:56:43 - Running 'TAPAS PhysIO Toolbox'
 
SPM12: spm_check_registration (v7759)              18:56:43 - 05/04/2022
========================================================================
Display /projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc2sub-39_T1w.nii,1
(all)  /projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc3sub-39_T1w.nii,1
 
Warning: fMRI volume and noise ROI have different orientation :
/projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/func/wasub-39_task-switching_run-04_bold.nii
 
/projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc2sub-39_T1w.nii
> In tapas_physio_create_noise_rois_regressors (line 98)
  In tapas_physio_main_create_regressors (line 338)
  In tapas_physio_cfg_matlabbatch>run_physio (line 1596)
  In cfg_run_cm (line 29)
  In cfg_util>local_runcj (line 1717)
  In cfg_util (line 972)
  In spm_jobman>fill_run_job (line 469)
  In spm_jobman (line 247)
  In preprocessing_task_data (line 54)

and then another one during coregistration:

05-Apr-2022 18:57:11 - Done    'Coregister: Estimate & Reslice'
05-Apr-2022 18:57:11 - Done
SPM12: spm_check_registration (v7759)              18:57:11 - 05/04/2022
========================================================================
Display /projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc2sub-39_T1w.nii,1
(all)  /projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc3sub-39_T1w.nii,1
 
Warning: Columns of X are linearly dependent to within machine precision.
Using only the first 526 components to compute TSQUARED.
> In pca>localTSquared (line 511)
  In pca (line 357)
  In tapas_physio_create_noise_rois_regressors (line 216)
  In tapas_physio_main_create_regressors (line 338)
  In tapas_physio_cfg_matlabbatch>run_physio (line 1596)
  In cfg_run_cm (line 29)
  In cfg_util>local_runcj (line 1717)
  In cfg_util (line 972)
  In spm_jobman>fill_run_job (line 469)
  In spm_jobman (line 247)
  In preprocessing_task_data (line 54)
Warning: fMRI volume and noise ROI have different orientation :
/projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/func/wasub-39_task-switching_run-04_bold.nii
 
/projects/chbh00080/data/derivatives/preprocessing/sub-39/ses-04/anat/wc3sub-39_T1w.nii
> In tapas_physio_create_noise_rois_regressors (line 98)
  In tapas_physio_main_create_regressors (line 338)
  In tapas_physio_cfg_matlabbatch>run_physio (line 1596)
  In cfg_run_cm (line 29)
  In cfg_util>local_runcj (line 1717)
  In cfg_util (line 972)
  In spm_jobman>fill_run_job (line 469)
  In spm_jobman (line 247)
  In preprocessing_task_data (line 54)

However when I check for the images orientation using SPM Check Reg and FLSeyes, the two images mentioned in the warning seem fine to me (see figures below)...any idea of what might cause the warning or what else I should check?

Thank you very much for your help!

Best wishes, Sara

[mrikasper]

Dear Sara,

It's great you are making continued use of the PhysIO Toolbox. I think the warning message about the orientation is not worrisome, you did the right thing, checking the registration in SPM (not sure whether FSL respects the header in the same way. I remember that in older releases of FSLview, they presented the image matrix, but not the transformations implied by the affine header).

Basically, PhysIO compares these affine headers in both nifti files and complains even about tiny differences, i.e..,

abs(Vroi.mat - Vimg(1).mat) > 1e-5)

So differences of about 10^-5 in any of the translation/rotation/shear etc. would trigger warning here, but might very well be indiscernable visually. What PhysIO then does is to coregister the images again with SPM, just to be sure. So if the final PhysIO plots overlaying the ROI and the functional image look reasonable to you, then all is good...unless you mind PhysIO performing this extra co-registration step. In a typical SPM pipeline, after warping to MNI space, the headers should all be equivalent though (because there has been reslicing). Did you use different pipelines to process the functional and structural (masking/tissue probability map) data here?

What I found more concerning is the other warning message you posted

Warning: Columns of X are linearly dependent to within machine precision.
Using only the first 526 components to compute TSQUARED.
> In pca>localTSquared (line 511)
  In pca (line 357)
  In tapas_physio_create_noise_rois_regressors (line 216)

This looks like there are identical images in your time-series (or duplicated voxels)...Did you do some motion censoring and scrubbing/interpolation of volumes, or maybe some global signal regression on this data? It could also be something harmless, e.g., all the background voxels were masked out and are zero time series. To be sure, you should carefully check that the output principal component time series (plot them via verbose.level =2) and the images of the PC weights (created files pc{01,02,etc}_scores_....nii) look plausible.

All the best, Lars

[SaraCal]

Dear Lars,

thank you very much for your prompt reply and the explanation. No I am using only one pipeline for the preprocessing (I am attaching it here so that you have all the info). Is there maybe something I need to change in the segmentation options?

preprocessing_task_data_job.m.zip

Regarding the second warning: no I don't do anything else with the data prior to the preprocessing... In the physio toolbox I am including the option to do censoring on volumes with excessive movement using MAXVAL, is that what might cause the warning?

I am attaching the 2 plots (for wm and csf) of the PCA components for one session and 2 screenshots from the images (wm and csf of PC 1 only - there is not much in the CSF ROI but they both look plausible to me?)

Thank you very much again!

Best wishes, Sara

[SaraCal]

Dear Lars, quick update: due to some messy Matlab paths I believe I was using the old version of PhysIO (from 2019) so I now changed to the latest version from 2021. Unfortunately I now get the following error:

Failed  'TAPAS PhysIO Toolbox'
Reference to non-existent field 'cm'.
In file "/projects/pbic1036/tapas-master/PhysIO/code/model/tapas_physio_create_noise_rois_regressors.m" (???), function "tapas_physio_create_noise_rois_regressors" at line 186.
In file "/projects/pbic1036/tapas-master/PhysIO/code/tapas_physio_main_create_regressors.m" (???), function "tapas_physio_main_create_regressors" at line 353.
In file "/projects/pbic1036/tapas-master/PhysIO/code/tapas_physio_cfg_matlabbatch.m" (???), function "run_physio" at line 1661.

Have you ever encountered this before? I am using the latest SPM12 version (version 7771)

Thank you again!

Best wishes, Sara

[SaraCal]

Dear Lars,

apologies for sending many messages, I have another update: I saw in another thread here that the issue with the non-existent field 'cm' can be solved by setting the verbose level to 1, so that is fixed now! Going back to the original issue (the warnings about image orientation), this is now the situation: the new version still gives me the warning below

Warning: fMRI volume and noise ROI have different orientation :
/projects/chbh00080/data/derivatives/preprocessing/sub-38/ses-01/func/wasub-38_task-switching_run-01_bold.nii
 
/projects/chbh00080/data/derivatives/preprocessing/sub-38/ses-01/anat/wc3sub-38_T1w.nii
> In tapas_physio_create_noise_rois_regressors (line 106)

and then proceeds to perform the coregistration (estimate & reslice) again. This extra step apparently creates another issue, because after the (extra) coregistration step there are very few voxels remaining in the functional file, for which then PhysIO issues the following error:

Error using tapas_physio_pca (line 33)
nVoxels <= nVolumes

I tried to debug this and I specifically noticed that before the extra coregistration step the number of voxels is fine (much greater than the number of volumes), but after that is drops to something like 46 voxels per volume...any idea of why that happens? I believe the extra coregistration step is done twice, once with the wc2 image and once with the wc3 image, and after the second time (coregistration to wc3) there are few voxels left. The two files that trigger the warning (as well as my preprocessing pipeline) are available in this drive folder in case it helps! https://drive.google.com/drive/folders/1owof9Xw2BP5UG18XW7am43C1-x0XVtjQ?usp=sharing

Best wishes, Sara

[SaraCal]

It appears as if there are few voxels in the noise ROI from the CSF? Should I change my noise modelling settings to see if I get more voxels there?

[benoitberanger]

Error using tapas_physio_pca (line 33) nVoxels <= nVolumes

This error pops when you have too few voxels in the mask, or long fmri runs (too many volumes), or both. It occurs to me especially with young subjects with small ventricles and a long task that cannot be split into several runs.

As you suggested, the solution is to be a bit loosen on the threshold for this mask, so more voxels are included. I would still keep at least an erosion of 1 voxel to ensure no partial volume effect. Please check the final mask after threshold and erosion, which should be written next the input mask.

[SaraCal]

Dear Lars, thanks a lot for your reply. I have now tested a few options and it all runs smoothly if I use 3 PCs and thresholds of 0.99 for WM and 0.85 for CSF (always cropping 1 voxel). Is 0.85 still a reasonable threshold for the CSF? :)

[mrikasper]

Dear Sara,

It was actually Benoît helping you out here - thank you, @benoitberanger! - because he also introduced the improved version on the noise ROIs extraction.

I agree with everything said, also the 85% threshold with 1 voxel erosion seems reasonable (but please check the mask image, as Benoît suggested).

The PCs you sent earlier look also fine, there seems to be some sudden motion(?) spikes at 150, 310 and 460 scans, do you see that in the realignment parameters as well? Could you send an updated version of these PC-plots, please?

What is still a bit puzzling is the need for coregistration, but it makes sense that if one tissue probability map has the wrong header, the other one will have so as well (wc and w3), because they stem from the same structural scan. I will have a look at the data you attached.

All the best, Lars

[mrikasper]

Dear Sara,

I just got your job with the data to run, thank you for providing it. There are some things to unpack:

1. The reason for the warning is trivial: the functional data has a lower resolution than the structural data (from which the mask is generated)...of course, we have to reslice to the same (functional) resolution then before we can extract the data.

   • This is not worrisome at all...maybe we should drop the warning in this case, what do you think, @benoitberanger? And rather put an info that the reslicing is necessary because of different voxel sizes...

2. I also got the warning of nVoxels < nVolumes - @benoitberanger , do you remember whether this was an issue before your reimplementation as well?

   • Conceptually, I would assume it should be possible to extract PCs also from rectangular matrices, if we have less voxels than volumes
   • Basically, one looks at the covariance of the signal over all voxels that one has. Maybe instead of the SVD of X with size [nVoxels, nVolumes] one has to perform a PCA on the covariance matrix (X'*X) (mean-centered) directly. There should still be enough samples (voxels) in that case. If you run

   [COEFF, SCORE, LATENT, EXPLAINED, MU] = tapas_physio_pca(Yroi, 'stats-pca');
   

   Matlab does offer a PCA even in the nVoxels < nVolumes case, and I think the toolbox-free version should be compatible with that. Our unit tests so far did not cover this case, that's why we missed it when integrating your implementation.

3. @SaraCal I have attached a version of tapas_physio_pca that switched to the Matlab Stats toolbox implementation in that case and also one of tapas_physio_create_noise_rois_regressors that avoids the error for verbose.level = 2 (which is better to set in order to get all the relevant output plots while tuning the noise rois parameters). gh162-noise-rois.zip

   • you have to overwrite the functions of the same name in the utils and model subfolder of spm12/toolbox/PhysIO.

4. @SaraCal, even after adjusting the threshold with 85% and 1 voxel cropping, the CSF mask does not contain a lot of voxels in the lateral ventricle

   • This is because the spm_erode cropping is very conservative, removing a voxel in all 3 dimensions, see here the difference in some transverse slices before/after cropping:
   • So having seen this, in your case it might be better to use a high threshold (like 0.99) and no cropping:
   • @benoitberanger , maybe we could crop in the (high-res structural) source space instead of the target space of the low-res functional images to have more control of how much one wants to crop

I hope that helps, Lars

[SaraCal]

Dear Lars and Benoît,

many thanks for your precious help (and apologies @benoitberanger for not addressing you earlier!) I will now implement the new functions you sent me. One final question: I guess the 0.99 threshold without cropping works best for this participant in particular, but should I just apply these parameters to all the subjects in my datasets or should I try to see if the others respond well to cropping 1 voxel?

Best wishes, Sara

[mrikasper]

Dear Sara,

I think having a consistent threshold and cropping over subjects would be beneficial, so you could try it for a few subjects (with verbose.level = 2 to visually assess the final (red) ROIs used for extraction. I think either a high threshold and now cropping or a more liberal threshold (even 80% is fine then) and 1 voxel cropping would be valid. The latter could have a higher partial voluming risk, while the latter reduces the statistics over voxels and might be more dominated by single voxel outliers.

All the best, Lars

[benoitberanger]

@SaraCal No worries.

@mrikasper

1. Warning at reslice Yes this warning can be transformed into a normal logging message, since it should not raise a particular attention.

2. nVoxels <= nVolumes This error message was introduced to ensure the user the PCA will be performed in the highly tested manner, which is with more channels than timepoints. This is a very conservative approach. Reading SPM VolumeOfInterest code (spm_run_voi + spm_regions), they indeed use the covariance matrix, and check the dimension of the input :

%-Compute regional response in terms of first eigenvariate
%--------------------------------------------------------------------------
if any(~isfinite(y(:)))
    error('Data contain NaN or Inf. Check the VOI definition.');
end
[m,n]   = size(y);
if m > n
    [v,s,v] = svd(y'*y);
    s       = diag(s);
    v       = v(:,1);
    u       = y*v/sqrt(s(1));
else
    [u,s,u] = svd(y*y');
    s       = diag(s);
    u       = u(:,1);
    v       = y'*u/sqrt(s(1));
end
d       = sign(sum(v));
u       = u*d;
v       = v*d;
Y       = u*sqrt(s(1)/n);

tapas_physio_pca could be adapted in a similar manner, after testing of course. Testing can be performed on synthetic input data.

4. Order for coregister+reslice, threshold, erode Swapping coregister+relsice and crop could make sense. This has to be tested first on real data. The development is minimal, but the testing is very demanding.

[benoitberanger]

Hello,

SPM team's strategy works to fetch one PC whatever the ratio nVoxel / nVolume, using the covariance matrix. I tested it on synthetic data. Here is almost the same code but that can fetch more than one PC :

% Each column is a voxel timeserie
[ nVolume , nVoxel ] = size(timeseries); % [ nVolume , nVoxel ]

not_finite = ~isfinite(timeseries);
if any(not_finite(:))
    warning('timeseries contains NaN or Inf, replacig it with 0')
    timeseries(not_finite) = 0;
end

% First regressor : mean timeserie of the ROI
mean_across_voxels = mean(timeseries,2); % [ nVolume , 1 ]

% Center data : remove temporal mean, mandatory step to perform SVD
mean_across_volumes = mean(timeseries,1); % [ 1 , nVoxel ]
timeseries_centered = timeseries - mean_across_volumes;

if nVolume > nVoxel
    [v,s,v] = svd(timeseries_centered' * timeseries_centered  );
    u       = timeseries_centered  * v/sqrt(s);
else
    [u,s,u] = svd(timeseries_centered  * timeseries_centered' );
    v       = timeseries_centered' * u/sqrt(s);
end

% Sign convention
d       = sign(sum(v)); 
u       = u.*d;
v       = v.*d;

% Singular values -> Eigen values
singular_values = diag(s);
eigen_values    = singular_values.^2/(nVoxel-1); % [ nVolume , 1 ]

% Eigen_values -> Variance explained
vairance_explained = 100*eigen_values/sum(eigen_values); % in percent(%) [ nVolume , 1 ]

% if nVolume >  nVoxel : [ nVolume , nVoxel  ]
% if nVolume <= nVoxel : [ nVolume , nVolume ]
PC = u*sqrt(s/nVolume); 

I can prepare a PR to solve a two issues:

• the warning at reslice
• the PCA when too few voxels, using the code above

Also I was thinking about adding a function that performs the threshold+erode and show it on a SPM on the fly, just like what is shown during the create_noise_regressor function. This could help the user to decide by visualizing dynamically the effect of their parameters.

About the order for coregister+reslice, threshold, erode, I didn't try it. However, after some thoughts, it think it will not be a good idea : the threshold+erode on the high-res mask will not eliminate enough voxels AFTER the reslice to fMRI grid.

Benoît