Several strategies have been proposed to model and remove physiological noise from resting-state fMRI (rs-fMRI) data, particularly at 7T, including contributions from respiratory volume (RV)¹ and heart rate (HR) fluctuations². Recent studies suggest that these fluctuations are highly variable across subjects and that their performance may benefit from an optimization at the dataset^3,4 or even voxel levels^2,5. Here, we systematically investigate the impact of the degree of specificity in the optimization of RV and HR models (group, dataset, region or voxel level), on the BOLD signal variance explained (VE) computed using a nested model approach, as well as the specificity and reproducibility of the functional connectivity (FC) measurements underlying three well-known resting-state networks (RSNs).

Data Acquisition

Two ~5mins scans were acquired in one session from 9 healthy subjects on a 7T whole-body scanner with a 32-channel receive RF coil. An SMS-EPI sequence was used with SMS-factor=3, TE=32ms, TR=2.5s, FA=75º, GRAPPA factor = 3, nominal echo spacing = 0.82ms, whole-brain coverage by 123 sagittal slices and 1.1mm isotropic resolution. Cardiac and respiratory data were simultaneously recorded using a pulse transducer (ADInstruments) and a pneumatic belt (UFI). A T1-weighted image was also acquired using multi-echo MPRAGE, with 1mm isotropic resolution⁶.

Data Pre-Processing

The rs-fMRI data pre-processing consisted of: slice timing correction, motion correction, slow drifts removal, periodic cardiac and respiratory fluctuations removal (using up to the 2^nd order RETROICOR) and spatial smoothing (3mm Gaussian kernel). The T1-weighted images were segmented into gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF), and co-registered with rs-fMRI and MNI images. Respiratory and cardiac signals were low-pass filtered and peak detection was performed on the filtered cardiac signal.

RV and HR models

The RV and HR contributions to the BOLD signal were modeled by shifting their time courses by a lag in the interval [-20 20]s in 1s steps, and then selecting either a Single-Lag (lag yielding maximum VE across GM) or a Dual-Lag (lags yielding the two first maxima VE across GM, one corresponding to the main positive peak and the other to the main negative peak) model⁷. Each model was optimized at four/three levels of specificity (for Single-Lag/Dual-Lag models, respectively), by selecting the optimal RV and HR time lags from the respective lag-VE curves as:

- Group: average across GM and subjects (1/2 optimal lags for the whole group)

- Dataset: average across GM in each dataset (1/2 optimal lags for each dataset)

- Region: average across three GM clusters in each dataset (3/6 optimal lags for each dataset), obtained by k-means clustering, within GM with k=3, on the voxelwise lag-VE curves (Fig.1)

- Voxel: the optimal lag for each voxel is used (because it is not possible to clearly identify a second VE maximum at the voxel level, Dual-Lag models are not used in this case)

RSN Functional Connectivity

A general linear model analysis was performed on the pre-processed and physiological noise-corrected rs-fMRI data using the average timecourses of three seed ROIs (PCC, IPS and SMA) orthogonalized to each other as regressors, in order to identify three well-established RSNs⁸. RSN group maps were obtained through higher-level analysis using a mixed-effects model and cluster thresholding (voxel Z>2.3 and cluster P<0.05). Individual FC measurements for each RSN were obtained as the respective GLM parameter estimates within the respective group map. FC specificity was evaluated by computing the ratio between the average FC inside each RSN map and the average FC across the whole brain. The voxelwise intraclass correlation coefficient (ICC) was computed for the FC maps to assess FC reproducibility.

Average VE results (Fig.2) show that more specific models are able to identify a higher putative spurious fraction of variance in the data, with dual-lag models outperforming single-lag models for group, dataset and region levels (repeated measures ANOVA, p<0.05); the former were therefore chosen for subsequent analysis in these cases. The most specific networks were obtained at the region level (Fig.3, Left), which exhibited the lowest reproducibility (Fig.3, Right). These results are in line with previous findings⁵, in that physiological noise fluctuations are highly reproducible within subjects, thus removing them results in lower ICC scores.Our results indicate that a region-based optimization of RV and HR physiological noise contributions is the most effective in removing spurious fluctuations from rs-fMRI data at 7T. Although voxelwise models explain more signal variance, they yield less specific FC measurements in RSNs, suggesting that they incur in overfitting to local fluctuations with no physiological meaning.