Introduction

fMRI has vast potential for insight into neural processes of the human spinal cord but is highly challenging in part due to the cord’s small size and proximity to physiological noise sources. Methods for physiological noise removal are well-developed in the brain but need further adaptation for the spinal cord. CompCor, a component based denoising technique for brain fMRI, defines a noise region-of-interest (ROI) in which principal component analysis (PCA) is used to identify nuisance regressors¹. Another method uses PCA in brain edge voxels to define nuisance regressors and improve removal of motion confounds in resting-state fMRI².

In an implementation of component-based denoising for resting-state spinal cord fMRI, a noise ROI is defined as anywhere that is “not-cord”; then PCA is used to identify nuisance regressors³. Components included are based on percent variance explained in the ROI, and noise is regressed out prior to motion correction to improve its efficacy. Here, we formalize and extend this method³ (herein called SpinalCompCor), examining the order of denoising relative to motion correction (Pipeline Comparison) as well as which principal components are included in denoising (Model Comparison).

Methods

Data collection:
Spinal cord fMRI was collected from 10 healthy participants (24.9±4.8y, 7F, 3M) on a 3T Siemens Prisma scanner. Data were acquired with GRE-EPI and ZOOMit selective excitation (TE/TR=30/2000ms, 1x1x3mm³, 25 slices, FA=90°, FOV=128x44mm²). Breath-hold (205 volumes) and hand-grip (300 volumes) task fMRI scans were each collected in two sessions (40 scans total).

SpinalCompCor PCA:
A mask was drawn around the spinal cord and cerebrospinal fluid (CSF) on the temporal mean fMRI volume, then dilated by 18 voxels. The original mask was subtracted from the dilation, to create the noise ROI (Fig. 1). Slice-wise PCA by eigenvalue decomposition was performed in the noise ROI. Following the automatic selection of PCs in the original implementation, the combination of PCs that represented up to 80% of the cumulative variance, or where the difference in variance between successive components was less than 5%, was calculated for each slice³. The mode number of PCs across slices meeting the criteria were designated for the Auto-Selected PCs model. The regressors were then either included in the final general linear model (GLM) or regressed from the data immediately.

Pipeline Comparison: Data were motion corrected by 2D slice-wise realignment⁴. SpinalCompCor was used to identify PC regressors. Data were registered to the PAM50 atlas⁵ and spatially smoothed⁶. In Pipeline A, SpinalCompCor identified PC regressors after motion correction and these regressors were included in the final GLM (Fig. 2). In Pipeline B, PC regressors were regressed out of the fMRI data first, followed by motion correction, and the final GLM did not include PC regressors. To compare, temporal signal-to-noise ratio (tSNR) of the corrected fMRI data (timeseries minus fitted variance) was calculated for each voxel and averaged across the spinal cord.

Model Comparison: Pipeline A was used to test variation in which PC regressors were included in the GLM. Ten expanded models were defined, voxel-wise F-Tests compared these to the reduced model (Fig. 3), and percent of significant voxels in the spinal cord was calculated for each dataset. The partial R² of the expanded models compared to the reduced model were calculated voxel-wise and averaged.

Results

Pipeline Comparison: The tSNR in the fMRI timeseries, after fitted variance was removed, was greater in Pipeline A compared to Pipeline B in 85% of datasets (Fig. 2).

Model Comparison: Voxels where the extended model fit significant additional variance (F-test, p<0.05, uncorrected) were distributed throughout the spinal cord/CSF region (Fig. 3). The percent of spinal cord voxels that were significant in models containing multiple PC regressors was significantly higher than with the original Auto-Selected PCs model, and with most models containing only one PC (Fig. 4A). Similar trends were found for the partial R² of each of the expanded models compared to the reduced (Fig. 4B).

Discussion

Applying SpinalCompCor after motion correction, versus before motion correction, improves tSNR in most datasets. This finding in 3T task-fMRI data is different from the original work in resting-state spinal cord fMRI at 7T³, in which component-based denoising was intended to improve the efficacy of motion correction algorithms. However, these implementations cannot be directly compared due to differences in data and protocol design.

Models including more PCs raise the percent significant voxels and partial R² values compared to the Auto-Selected PCs. Both metrics indicate that including the second PC creates a better model. As the included number of PC regressors increases, the F-Test percent of significant voxels begins to plateau. These observations suggest that multiple PC regressors are warranted in the GLM, but there is a point of diminishing returns.

This method is based on the central assumption that signal fluctuation in the “not-cord” ROI can be attributed to motion confounds or physiological noise and does not represent signals originating from neurovascular coupling. Future work will focus on the origin of signal fluctuation in the “not-cord” ROI, task-correlated PCs, and spectral analyses.

Conclusion

SpinalCompCor uses PCA to derive nuisance regressors from a region outside of the spinal cord, and these regressors contribute to an improved denoising model when included in spinal cord fMRI analysis.

Acknowledgements

Research supported by the Craig H. Neilsen Foundation (595499). K.J.H. is supported by an NIH funded Training Program (T32EB025766).