Towards reliable spinal cord fMRI: Assessment of common imaging protocols

Functional magnetic resonance imaging (fMRI) has revolutionized the investigation of brain function. Similar approaches can be translated to probe spinal mechanisms. However, imaging the spinal cord remains challenging, notably due to its size and location. Technological advances are gradually tackling these issues, though there is yet no consensus on optimal acquisition protocols. In this study, we assessed the performance of three sequences during a simple motor task and at rest, in 15 healthy humans. Building upon recent literature, we selected three imaging protocols: a sequence integrating outer volume suppression (OVS) and two sequences implementing inner field-of-view imaging (ZOOMit) with different spatial and temporal resolutions. Images acquired using the OVS sequence appeared more prone to breathing-induced signal fluctuations, though they exhibited a higher temporal signal-to-noise ratio than ZOOMit sequences. Conversely, the spatial signal-to-noise ratio was higher for the two ZOOMit schemes. In spite of these differences in signal properties, all sequences yielded comparable performance in detecting group-level task-related activity, observed in the expected spinal levels. Nevertheless, our results suggest a superior sensitivity and robustness of patterns imaged using the OVS acquisition scheme. To analyze the data acquired at rest, we deployed a dynamic functional connectivity framework, SpiCiCAP, and we evaluated the ability of the three acquisition schemes to disentangle intrinsic spinal signals. We demonstrated that meaningful subdivisions of the spinal cord's functional architecture could be uncovered for all three sequences, with similar spatio-temporal properties across acquisition parameters. Cleaner and more stable components were, however, obtained using ZOOMit sequences. This study emphasizes the potential of fMRI as a robust tool to image spinal activity in vivo and it highlights specificities and similarities of three acquisition methods. This represents a key step towards the establishment of standardized spinal cord fMRI protocols.


Introduction
The spinal cord is a pivotal part of the central nervous system (CNS), at the interface between the brain and the periphery. As such, it is critical to support sensorimotor functions, and it can be involved in many neurological disorders (e.g., multiple sclerosis , spinal cord injuries ( Silva et al., 2014 ), neuropathic pain ( Colloca et al., 2017 ) or amyotrophic lateral sclerosis ( Hardiman et al., 2017 ). Despite this primary importance, our view of the CNS has for long been mostly focused on the brain, while the spinal cord has so far received limited attention. This dearth of studies partly pertains to the limited availability of methods to efficiently probe spinal cord function in humans . Although knowledge on that topic can be gained from animal models, spinal circuits are not entirely preserved or temporal resolution -as well as on limiting artifacts. In particular, movements of external tissues (e.g., chest, jaw, etc.) can be detrimental as they introduce temporal variations in field inhomogeneity, and thus lead to distortions and blurring of the images along the phase-encoding (PE) direction ( Saritas et al., 2014 ).
Since the seminal work of Yoshizawa in 1996( Yoshizawa et al., 1996, a number of studies have implemented different methodologies to circumvent these impediments (see Eippert et al., 2017a;Powers et al., 2018;Stroman et al., 2014;Summers et al., 2014;Wheeler-Kingshott et al., 2014 for reviews). Nevertheless, the heterogeneity of the employed strategies has, so far, precluded the emergence of standardized guidelines. Notably, most of the early studies were conducted using non-standard acquisition techniques, favoring a contrast mechanism termed signal enhancement by extravascular water protons (SEEP) ( Stroman et al., 2001a( Stroman et al., , 2001b, in place of the well-established BOLD contrast. The SEEP technique, relying on T2-weighted spin-echo MRI, was suggested to provide a better image quality. The findings related to this contrast mechanism are, however, controversial and several groups unsuccessfully tried to reproduce these results ( Bouwman et al., 2008 ;Jochimsen et al., 2005 ;Moffitt et al., 2005 ). Nowadays, T2 * -weighted gradient-echo-planar imaging (EPI), prevalent in brain fMRI to image BOLD signal, is most commonly used for spinal cord fMRI, though approaches such as multi-shot 3D fast-field echo imaging have recently shown promising results ( Barry et al., 2021 ).
Although gradient-echo EPI allows rapid imaging and high signal-tonoise ratio, it is particularly sensitive to susceptibility variations, prominent in the spinal cord due to the various surrounding tissues (Powers et al., 2018) . Several techniques have been proposed to mitigate these artifacts and, thus, to maximize signal quality. For instance, spinal cord imaging can benefit from methods limiting the extent of the field-ofview (FOV) in the PE direction (i.e., anterior-posterior), so that signals from external tissues are not included. A simple implementation for this is to use a standard sequence with outer volume suppression (OVS), which consists in the application of saturation bands located externally to the cord ( Wilm et al., 2007 ). This technique was employed in a number of fMRI studies, with promising results ( Eippert et al., 2017b( Eippert et al., , 2009Geuter and Büchel, 2013 ;Kong et al., 2014 ;Sprenger et al., 2012 ;van de Sand et al., 2015 ). Alternatively, inner FOV imaging can be achieved through dynamic pulses specifically exciting the region-of-interest in the PE direction. This method is now implemented in most scanners under different denominations (e.g., ZOOMit for Siemens, FOCUS for GE or iZOOM for Philips). Since encoding is limited to the targeted region, it allows for reduced acquisition time and high spatial resolution, while reducing susceptibility artifacts. As a matter of fact, the applicability of inner FOV approaches has been widely demonstrated for diffusion weighted imaging of the spinal cord ( Alizadeh et al., 2017;Finsterbusch, 2012;Samson et al., 2016;Saritas et al., 2008;Wheeler-Kingshott et al., 2002;Zaharchuk et al., 2011 ). More recently, several fMRI studies also underlined their potential to image spinal cord function ( Kinany et al., 2020( Kinany et al., , 2019Weber et al., 2020Weber et al., , 2016aWeber et al., , 2016b. Notwithstanding the compelling results, both using OVS and inner FOV techniques, there is yet no consensus as regards optimal acquisition parameters, generalizability across sequences, nor differences in performance. Addressing these questions can, therefore, offer a stepping stone to larger multicentric studies aiming to establish a reproducible pipeline to image spinal cord function , similarly to recent efforts in proposing an accessible and reliable pipeline for quantitative MRI to probe spinal cord structure ( Cohen-Adad et al., 2021a, 2021b. To this end, we systematically compared three common spinal cord fMRI protocols during a simple motor paradigm (wrist adduction), previously demonstrated to elicit localized activity ( Kinany et al., 2019 ), and at rest. Specifically, we selected one standard sequence that uses outer volume suppression (OVS), as well as two inner FOV implementations (ZOOMit), with different spatiotemporal resolutions. Capitalizing on state-of-the-art processing tools tailored to the spinal cord (e.g., physiological noise correction, common template, slice-wise motion correc-tion, etc.) ( De Leener et al., 2017 ;Eippert et al., 2017 a), we evaluated the suitability of the three sequences to detect task-evoked and spontaneous activity, using analytical approaches on par with brain standards. Our results allowed us to assess the characteristics of these acquisition schemes, while probing the robustness of spinal cord fMRI results against variations in the acquisition pipeline. These investigations represent an important step towards the development of a reliable and generalizable pipeline for spinal cord fMRI.

Participants
Sixteen healthy volunteers were enrolled in this study. All subjects had normal or corrected-to-normal vision and no history of neurological or motor disorders. One participant had to be excluded due to technical problems (misalignment of the imaging stack). The final sample consisted of 15 participants (8 female, 28.9 ± 4.2 years old). The study had been approved by the Commission Cantonale d'Éthique de la Recherche Genève (CCER, Geneva, Switzerland, 2016-01,566). All participants gave their written informed consent to participate.

Experimental paradigms
During the experiment, each participant performed two task runs and one rest run for each acquisition scheme. During the task runs, subjects were asked to execute a bilateral wrist adduction movement (i.e., wrist towards the outside, with an ulnar deviation). In a previous study, we had demonstrated through a combination of electromyographic and fMRI recordings that this movement elicited localized spinal activity (Kinany et al., 2019) . Within each run, movements were performed in blocks of 18 s (8 blocks of rest alternated with 8 blocks of movement and a final block of rest) and an entire run lasted approximately five minutes. Instructions were displayed on a screen (fixation cross ' + ' during the rest blocks and text indicating 'movement repetitions' during the task blocks). Auditory cues were used to inform the subject of the different phases, as well as to ensure a uniform rhythm across participants (10 movements per block). Rest runs (i.e., no overt task) lasted 7 min and 30 s. For all subjects, the experiment was carried out with the following structure: Rest A / Run 1A -Run 1B -Run 1C / Rest B / Run 2A -Run 2B -Run 2C / Rest C. A, B, C corresponded to the three types of acquisition and were randomized over subjects.

Data acquisition
All experiments were performed on a Siemens 3T Prisma scanner (Erlangen, Germany), equipped with a 64-channel head (only inferior element, HC7, was used) and neck coil (both anterior and posterior elements, NC1 and NC2, were used -i.e., 24 channels), as well as a spine coil (only the upper element of the spine coil, SP1, was used, following the optimal coil combination defined by the scanner). Subjects were placed in the scanner in supine position, using cushions and a soft cervical collar to stabilize the neck and minimize its curvature.
Functional images were acquired using T2 * -weighted gradient echoplanar imaging (EPI) sequences, using different sets of acquisition parameters. Based on the current spinal cord fMRI literature, we opted for three acquisition schemes: one standard sequence integrating OVS (hereafter: OVS), based on Eippert et al. (2017b) , and two ZOOMit inner field-of-view sequences with different spatial and temporal resolutions (hereafter: Z3mm and Z5mm), similar to Kinany et al. (2019) . These protocols relied on built-in Siemens sequences. Details of the parameters are provided in Table 1 . For OVS and Z3mm, parameters were kept identical to the reference sequences (see Eippert et al. 2017b, Kinany et al. 2019. Z5mm, on the other hand, was adapted from Z3mm. In order to minimize signal dropout , we thus decided to favor the use of a shorter TE, which was not possible for Z3mm using our scanner. OVS was performed thanks to spatially-selective saturation pulses that were applied above and below the imaging stack to reduce sensitivity to flow effects and flow rephasing in slice direction (Fig. S1A). Additional saturation pulses anterior and posterior to the target region were also used to reduce pulsatile blood flow artifacts. Imaging stacks were positioned to cover a region spanning the C5 to C8 spinal levels (Fig. S1B). For all functional acquisitions, slices were positioned perpendicularly to the spinal cord and aligned to the intervertebral disks (Fig. S1C). As EPI is prone to signal losses and distortions due to magnetic field shifts (Powers et al., 2018) , efficient shimming of the magnetic field is crucial to limit such artifacts. For this reason, the shim volume was manually set to be focused on the spinal cord ( Eippert et al., 2017a;Ellingson and Cohen-Adad, 2014 ) (Fig. S1A). Before undergoing processing, all functional images were visually inspected and bottom slices with insufficient signal were removed (between 0 and 5 slices over all subjects and runs).
A high-resolution T2-weighted anatomical image was also acquired in all participants using a SPACE sequence (single slab 3D turbo spin echo sequence with a slab selective, variable excitation pulse, TR = 1500 ms, TE = 135 ms, echo train length = 74, flip angle = 140, resolution = 0.4 × 0.4 × 0.8 mm, sagittal orientation). The imaged region extended from the upper cervical to the upper thoracic spine.

Motion correction
To ensure proper volume alignment within and between runs, a two-step motion correction procedure was employed: i) For each run, slice-wise realignment was performed using the mean functional image as a reference, with spline interpolation and a mean-square cost function, using the Spinal Cord Toolbox (SCT) ( De Leener et al., 2017 ). Computations were limited to a cylindrical mask automatically drawn along the spinal cord, to exclude areas outside the vertebral column that may move independently. ii) All runs were then aligned to the rest run with FLIRT (FMRIB's Linear Image Registration Tool), using six degreesof-freedom rigid-body transformations with spline interpolation and a normalized correlation cost function ( Jenkinson et al., 2002 ). Motion parameters from step i) were used to compute the mean (i.e., average over slices and volumes) framewise displacement (FD). While these steps can reduce the noise in BOLD time series, severe artifacts due to large volume-to-volume movements may persist ( Caballero-Gaudes and Reynolds, 2017 ;Parkes et al., 2018 ). In order to account for this, motion scrubbing was also employed ( Power et al., 2014 ). Noise regressors corresponding to the artifactual frames were obtained based on variations in image intensity between volumes, using FSL's outlier detection tool (DVARS metrics: root mean square intensity difference of volume N to volume N + 1, within the spinal cord mask) with a box-plot cutoff (75th percentile + 1.5 x the interquartile range).

Segmentation
Using the SCT's function sct_deepseg_sc ( Gros et al., 2019 ), a spinal cord binary segmentation mask was automatically drawn based on T2weighted anatomical images. As for functional images, the spinal cord was manually segmented, for each sequence, using the mean motioncorrected functional images (average of the three aligned runs). Similarly, a binary mask of the cerebrospinal fluid (CSF) was manually drawn, keeping a gap of one voxel with the spinal cord. In order to quantify the contrast between the spinal cord and CSF in functional images, we assessed the absolute difference between the mean signals in the white matter (to exclude influence from non-neighboring gray matter voxels) and the CSF (removing the 20th lower percentile of voxel intensities, to account for the presence of nerves roots), divided by the standard deviation of the white matter signal .

Physiological noise modeling
In the spinal cord, physiological noise can be particularly pronounced, due to the surrounding respiratory and cardiac organs ( Piché et al., 2009 ;Vannesjo et al., 2018 ). To limit the impact of these confounds, physiological recordings were performed using a photoplethysmograph and a respiratory belt (Biopac MP150 system, California, USA), synchronized with the scanner triggers. Using these signals, we adopted an approach based on the RETROspective Image CORrection (RETROICOR) procedure ( Glover et al., 2000 ) to model physiological noise. Briefly, this method assumes physiological processes to be quasiperiodic, such that cardiac and respiratory phases can be uniquely assigned for each volume. Signals can then be modeled using a low-order Fourier expansion based on the phases at the time of image acquisition. Following recommendations for spinal cord fMRI, we included higher order Fourier terms, as well as terms accounting for the interactions between cardiac and respiratory cycles (Brooks et al., 2008 ;Kong et al., 2012 ). Specifically, 32 slice-specific noise regressors were extracted, for each run, using the FSL's physiological noise modeling (PNM) tool. To account for the noise present in the CSF, an additional regressor was computed using the CSF signal (slice-wise averaged time courses from the 10% of CSF voxels showing the highest variance).

Spatial normalization
Functional images were normalized to the PAM50 template ( De Leener et al., 2018 ) following a two-step transformation procedure using the SCT (De Leener et al., 2017) : (i) Anatomical-to-template: landmarks for each vertebra were manually generated and used, together with a binary mask of the spinal cord (see 2.4.2), to perform nonrigid registration to the PAM50 template. (ii) Functional-to-anatomical: functional volumes were registered to the corresponding anatomical scan, using binary masks of the spinal cord (see 2.4.2 ). Finally, the warping fields obtained in both steps were concatenated to obtain the functional-to-template transformation. Mean normalized anatomical and functional images are presented in Fig. S2.

Data analysis
The different steps of the analysis performed after data preprocessing are shown in Fig. 1 . In a first step (see 2.5.1 ), we investigated the signal quality for the three sequences. We then probed the sequences' ability to reliably detect task-related activation (see 2.5.2 ). Finally, we characterize their potential to capture the intrinsic functional architecture of the spinal cord at rest (see 2.5.3 ).

Signal-to-noise ratio
To inform on signal quality for each acquisition scheme, the entire rest runs (i.e., 7 ′ 30'') were used to evaluate the temporal (tSNR) and spatial (sSNR) signal-to-noise ratio  . Those two measures allow to assess the data quality over time or over voxels, respectively. Specifically, the tSNR was obtained from the motion corrected time-series. For each voxel, we computed the average signal intensity across time, divided by its standard deviation. To allow comparison between sequences, tSNR values of each subject were averaged within the spinal cord mask (C5 to C8) and compared using repeated measures ANOVA. On the other hand, the sSNR was computed on the mean motion corrected image (i.e., mean over the 7 ′ 30''), as the average signal intensity across voxels, in the spinal cord mask, divided by its standard deviation. Again, values were compared between sequences using repeated measures ANOVA. To further explore signal quality, mean tSNR and sSNR values in gray matter and white matter masks were also extracted independently. This computation was performed in the native space. To this end, masks were transformed from the PAM50 template space to each individual space and binarized with a threshold of 0.2. Finally, for each sequence, paired t-tests were used to compare values derived from the two tissue types, for each acquisition scheme (Bonferroni corrected for multiple comparison, with n = 6).

Task-related analysis
Task-related analysis was carried out using a general linear model approach (GLM), implemented in FSL's fMRI Expert Analysis Tool (FEAT). At the subject-level, the preprocessed images of each run (in the native space, after motion correction) were first denoised using FSL's fMRI Expert Analysis Tool (FEAT) to regress out the physiological noise regressors, the motion parameters, as well as the outlier volumes (see 2.4.1 ). Denoised time series were then normalized to the PAM50 template (De Leener et al., 2018) and spatially smoothed using a 3D Gaussian kernel with a full width half maximum (FWHM) of 2 × 2 × 6 mm 3 . High-pass temporal filtering (cut-off frequency = 90 s) was applied on the obtained time-series, for each task run, which were fed to a firstlevel statistical analysis using FMRIB's Improved Linear Model (FILM) with local autocorrelation correction ( Woolrich et al., 2001 ). The design matrix included the explanatory variables, modeled based on the temporal dynamics of the task blocks and convolved with three optimal basis functions using FMRIB's Linear Optimal Basis Set (FLOBS) ( Woolrich et al., 2004b ), in order to account for differences of hemodynamic response function (HRF) between regions and subjects. The second and third waveforms (i.e., the temporal and dispersion derivatives, respectively) were orthogonalized to the first waveform. For each subject and movement, the parameter estimates corresponding to the first FLOBS waveform (i.e., task against rest), obtained independently for the two runs, were combined by passing them up to a second-level analysis, in which task-specific (i.e., for each acquisition scheme) subject level activation maps were obtained using a fixed-effects model.
Parameter estimates from these second-level analyses were finally fed to a third-level analysis aimed to yield group activation maps, using FMRIB's Local Analysis of Mixed Effects (FLAME) stages 1 and 2 with outlier detection ( Woolrich, 2008 ;Woolrich et al., 2004 a). Z statistic images were thresholded using clusters determined by Z > 2.3 and a cluster-defining threshold of p < 0.001 to account for multiple comparisons, in line with current recommendations ( Eklund et al., 2016 ). Group-level analyses were performed within a common mask, spanning C5 to C8 spinal levels, for the three acquisition schemes. The localization of the activity was probed by counting the number of voxels (in the PAM50 template space) in each spinal segment (C5 to C8), as well as within each hemicord: left, right, dorsal (i.e., posterior) and ventral (i.e., anterior).
For control purposes, preprocessed rest runs were temporally cropped to match the duration of task runs and analyzed following the same procedure as the task runs. Group average activation maps were obtained for each sequence (thresholded using clusters determined by Z > 2.3 and a cluster-defining threshold of p < 0.001). For each subject, acquisition scheme and run, the number of active voxels and their average Z-score were also computed. To this end, individual Z-statistic images were normalized to the PAM50 template and thresholded ( Z > 2.3, uncorrected). These control metrics were then compared to values obtained for the task runs of each sequence, using two-tailed paired ttests and Bonferroni correction to account for multiple comparisons. To enable comparison with individual rest runs, note that the task runs were considered independently (second-level statistics) and averaged for each subject.
Finally, we explored repeat reliability (i.e., between-run similarity). To this end, group-level activation maps were obtained independently for Runs 1 and Runs 2, for each sequence. Their similarity was then assessed by means of Pearson correlation coefficients. For comparison purposes, the similarity with the group-level activation maps of the corresponding rest runs was also computed.

Resting-state analysis
The resting-state scans were processed as presented in Section 2.4 and denoised (regression of physiological noise, motion parameters and outlier volumes) using FSL's fMRI Expert Analysis Tool (FEAT). The denoised time-series were then normalized to the PAM50 template (De Leener et al., 2018) and spatially smoothed using a 3D Gaussian kernel with a FWHM of 2 × 2 × 6 mm 3 .
Following the SpiCiCAP framework ( Karahano ğlu and Van De Ville, 2015 ; Kinany et al., 2020 ), we used the denoised time-series (native space) and applied the Total Activation framework (TA) ( Karahano ğlu et al., 2013 ) to obtain activity-inducing signals (i.e., the signal deconvolved from the hemodynamic blur). Innovation signals were computed as the temporal derivative of these activity-inducing time-series to represent transient activity. A two-step thresholding process was used to highlight significant innovations: (i) temporal thresholding: for each voxel, TA was applied on phase randomized data to yield a surrogate distribution, on which the 5% confidence interval was used to select significant voxels, (ii) spatial thresholding: only the innovation frames with more than 5% of active voxels were kept.
K-means clustering was performed for each sequence individually, using the significant innovations obtained from all subjects. For the sake of comparability across sequences, we selected the same number of iCAPs (K) for all sequences. In order to ensure that the selected levels of details enabled reliable partitioning of the data, we systematically evaluated the reproducibility of the clustering for different values of K (Fig.  S5). With a subsampling scheme, we performed clustering using different K in 100 random subsets of the data (10 subjects). K values were chosen to span a low-granularity range (from 4 to 20, in steps of 4, corresponding to the number of spinal levels) and a high-granularity range (from 20 to 60, in steps of 10). Each subsample clustering solution was compared to the global clustering obtained with the 15 subjects, with the adjusted mutual information (AMI) ( Vinh et al., 2009 ). This metric reflects the similarity of two discrete assignments (i.e., by comparing the assignments of the significant innovation frames to the different clusters) and is corrected for the effect of chance, thus limiting bias towards large numbers of clusters. Values range from 0 (chance level) to 1 (equal partitions). While different trends were observed for the three acquisition schemes, K = 4 (low-granularity) and K = 30 (high-granularity) appeared as reliable solutions and these numbers were selected for subsequent steps. For all non-noisy iCAPs, inter-subject stability was evaluated using the mean cosine similarity over all pairs of subject-specific maps (computed as the mean over the frames of this subject assigned to an iCAP) for a particular iCAP.
To investigate iCAPs' spatial properties, the thresholded iCAP maps ( Z > 1.6 for K = 4 and Z > 5 for K = 30) were used to compute, for each component, the proportion of voxels found in the different spinal levels or atlas regions (Lévy et al., 2015) . Using a hard assignment procedure, fine-grained iCAPs were then uniquely matched to the atlas region containing the maximum number of voxels. The accuracy of this assignment was confirmed using Dice coefficients ( Dice, 1945 ). In the case of bilateral iCAPs (i.e., less than 70% of their voxels in a single hemicord), the matching procedure was applied independently for the two hemicords.
Finally, we obtained subject-specific time courses by regional averaging of the activity-inducing signals for each iCAP. These time courses were Z-scored and thresholded (|Z| > 1) to identify strongly active and deactive time points. Temporal properties were summarized using the average and total durations of each iCAP. Temporal overlap between pairs of iCAPs was explored by reporting Jaccard indices (i.e., percent joint activation time), in order to compute couplings and anti-couplings. Fig. 2 A presents the mean functional images after motion correction from one representative subject. For all three acquisition schemes, the characteristic gray matter butterfly could be observed in most slices. Interestingly, the spinal cord to CSF contrast differed across sequences ( p < 0.001, repeated measure ANOVA). Images had a significantly stronger contrast ( p < 0.01, post-hoc t -test with Bonferroni correction) for Z3mm (contrast = 1.71 ± 0.22, mean ± standard error (SE)) than for the sequences with a larger slice thickness (contrast = 0.95 ± 0.15 for OVS and 0.73 ± 0.11 for Z5mm). Furthermore, OVS and Z5mm appeared more prone to signal losses, mostly occurring on the dorsal side of the spinal cord. A quantitative assessment of tSNR and sSNR further revealed different signal properties for the three sequences ( Fig. 2 B). tSNR was significantly higher ( p < 0.001) for the OVS acquisition scheme (10.06 ± 0.32, mean ± standard error (SE)) compared to Z3mm (6.80 ± 0.20) and Z5mm (7.45 ± 0.13). An opposite trend was observed for sSNR, with larger values ( p < 0.001) for Z3mm (6.10 ± 0.31) and Z5mm (5.33 ± 0.47), compared to OVS (3.25 ± 0.13). Of note, for all three sequences, tSNR was not uniform throughout the spinal cord, as periodic signal losses were present around intervertebral disks ( Fig. 2 C). Moreover, both sSNR and tSNR values were significantly higher in the gray matter compared to the white matter (Fig. S3). The mean frame displacement (FD) showed no significant difference within-sequence (i.e., for the different runs) (see Table S1). Nevertheless, OVS, Z3mm and Z5mm had, respectively, an average FD of 0.36 ± 0.04 mm (mean ± SE), 0.14 ± 0.04 mm and 0.23 ± 0.05 mm, highlighting distinct movement properties across sequences ( p < 0.001, repeated-measures ANOVA followed by paired t-tests, Bonferroni corrected for multiple comparison). In particular, examination of the temporal and spectral properties of the motion parameters suggested that the increased motion observed using OVS was due to fluctuations in the range of the respiratory rhythm ( ∼0.3 Hz, ( Fleming et al., 2011 ) (Fig. S4).

Group-level activation maps
Although image quality (i.e., sufficient contrast, limited distortion, etc.) and SNR are important factors to evaluate signal integrity, they do not fully inform on BOLD sensitivity. To evaluate this, we investigated task-related activity using the GLM approach. Significant task activity was detected for all three acquisition schemes ( Fig. 3 A, B). The extent of activity was smaller for Z5mm (1885 active voxels) compared to OVS and Z3mm (3191 and 3965 active voxels, respectively). In a previous study (Kinany et al., 2019) , we showed through electromyographicbased and fMRI-based spinal maps that activity related to wrist adduction is mainly focused in C7 and C8, which is expected from anatomical knowledge ( Kendall et al., 2010 ). Importantly, activation patterns imaged in the current study coincided with this region ( Fig. 3 B, C). In particular, 84.1% of active voxels fell into spinal levels C7 and C8 when using OVS, 80.2% for Z3mm and 86.7% for Z5mm. Finally, group-level activations for the three acquisition schemes were distributed over the different hemicords, dorso-ventrally and laterally, as expected for bilateral dynamic movements ( Fig. 3 C).  Group activation maps (mixed-effects modeling) of the control analysis (i.e., rest runs), for the three sequences. Maps are thresholded at a Z-score > 2.3 (cluster-defining threshold of p < 0.001) and normalized to the PAM50 template (De Leener et al., 2018) . Coronal views are presented. n.s. indicates that no significant activity was detected. B. Comparison of the number of active voxels (left panel), and mean Z-score (right panel) during the control and task runs. Each colored line corresponds to one subject. The black line corresponds to the mean ± SE. * corresponds to p < 0.01 and * * to p < 0.001 (two-tailed paired t -test, Bonferroni corrected for multiple comparisons). C = control, T = task.

Control analysis
To highlight the genuine nature of these activation patterns, we performed a control analysis in which rest runs were analyzed following the same pipeline. As expected, no voxel exceeded the significance threshold at the group-level for all three sequences ( Fig. 4 A). When probing the activity in each individual subject, we reported activations both in task and control runs. Nevertheless, the spatial extent of these activations, as illustrated by the number of active voxels, was significantly lower during the control runs (mean ± SE for OVS, Z3mm and Z5mm, respectively: 311 ± 60, 603.3 ± 93.6 and 565.1 ± 69.2) than during the task runs (2378.7 ± 492.9, p < 0.01, 988.6 ± 147.6, p < 0.05 and 1198.4 ± 180.7, p < 0.05) ( Fig. 4 B). The magnitude of activity (average Z score) was significantly lower during the control runs for OVS (2.59 ± 0.04 vs 2.83 ± 0.06, p < 0.05), but not for Z3mm (2.67 ± 0.03 vs 2.72 ± 0.03) and Z5mm (2.66 ± 0.03 vs 2.73 ± 0.03).

Repeat reliability
Taking advantage of the two runs acquired for each subject and sequence, we then evaluated the repeat reliability of the task-related activation patterns. Maps obtained with the OVS sequence appeared more stable (Pearson correlation = 0.47) than those obtained with Z3mm ( = 0.15) and Z5mm ( = 0.14). However, for the three sequences, these correlations were larger than the ones obtained between each task run and the corresponding rest run. The latter were, indeed, weakly negative ( = − 0.07 for OVS, − 0.4 for Z3mm and − 0.03 for Z5mm, average over the correlations between the two task runs and the rest run).

Spatial organization
We used the SpiCiCAP framework to extract iCAPs for two levels of granularity: K = 4 (low-granularity) and K = 30 (high-granularity).

Neuroanatomical relevance
Considering the localized organization of fine-grained iCAPs, we investigated their spatial properties by probing the distribution of their voxels in anatomically-defined atlas regions ( Lévy et al., 2015 ). ICAPs spatial properties appeared to be focal and in line with the underlying neuroanatomy. Indeed, iCAPs obtained using the three sequences could be matched to an atlas region with a high accuracy (Dice coefficients, mean over iCAPs ± SD, OVS: 0.54 ± 0.12, Z3mm: 0.57 ± 0.10 and Z5mm 0.57 ± 0.09) ( Fig. 5 B). There was no significant difference between the Dice coefficients of the different sequences ( p = 0.35, one-way ANOVA).
To further explore the neuroanatomical nature of fine-grained iCAPs, we evaluated the number of iCAPs assigned to each atlas region ( Fig. 6 ). When focusing on iCAPs with a unilateral cluster (dark colors in Fig. 6 ), the overall distribution was similar for the three acquisition schemes. As expected from our previous observations ( Kinany et al., 2020 ), a large portion of iCAPs were found along the cortico-spinal tract pathway (CST, 5 iCAPs for OVS, 10 for Z3mm and 9 for Z5mm) and the dorsalcolumn-medial lemniscus pathway (DCML, 7 iCAPs for OVS, 9 for Z3mm and 11 for Z5mm). In addition to regions associated with these major neural pathways, iCAPs were also corresponding to extrapyramidal descending pathways (vestibulospinal, rubrospinal and tectospinal tracts) (2 iCAPs for OVS, 4 for Z3mm and 8 for Z5mm). These pathways are formed by smaller tracts originating in the brainstem, and are involved in unconscious muscle control (balance, posture, etc.). ICAPs linked to the spino-olivary and spinocerebellar tracts, ascending routes conveying unconscious proprioceptive information to the cerebellum ( Rea, 2015 ), were also reported (6 iCAP for OVS). Finally, the remaining unilateral iCAPs were found in the ventral and intermediate parts of the gray matter (1 for OVS and 5 for Z3mm). As for bilateral ICAPs (light colors in Fig. 6 ), regions were similar to the ones observed in components with a single cluster. Specifically, iCAPs included bilateral clusters in the dorsal column (2 for OVS, 1 for Z3mm and 1 for Z5mm) as well as between motor regions (7 iCAPs containing CST and vestibulospinal regions for OVS, 1 iCAP containing CST and tectospinal regions for Z3mm, and 1 iCAP containing bilateral ventral horns for Z5mm).
To examine the nature of these overlaps, we divided them based on their signs (couplings or anti-couplings) and their locations (within or between spinal levels) ( Fig. 7 D). A three-way ANOVA revealed a significant main effect of the sign ( p < 0.001) and location ( p < 0.001), as previously suggested ( Kinany et al., 2020 ), but no main effect of the sequence ( p = 0.70). Further decomposing this 3-way interaction, we highlighted a significant two-way interaction between sign and location ( p < 0.001). Post-hoc t-tests (Bonferroni corrected for multiple comparisons) confirmed that couplings were stronger than anti-couplings ( Fig. 7 D), and most prominent within-compared to between levels. In contrast, no significant difference was reported between anti-couplings at different locations.

Discussion
Spinal cord fMRI is a promising approach to shed light on human spinal mechanisms in vivo. Yet, the progression of the field has been relatively limited, with technical hurdles hindering its deployment. As MRI-based technology advances (e.g., higher field strength, innovative pulse sequences or improved analyses), new tools are gradually tackling these limitations. Nevertheless, current approaches remain heterogeneous, and a systematic characterization of the methodologies is still lacking. Here, we conducted a comparative study to evaluate three commercially available acquisition protocols implementing pulse sequences based upon recent literature. We showed that, regardless of differences in signal properties, all three sequences could be used to efficiently image motor-evoked activity and spontaneous fluctuations.
The first phase of our analyses focused on evaluating image and signal quality for the three acquisition schemes. This emphasized their dis- Fig. 5. iCAPs spatial patterns. For the three sequences, spatials maps of the low-granularity iCAPs ( A. , coronal and axial views) and the high-granularity iCAPs ( B. , axial views, coronal views are presented in Fig. S6) are presented. Thresholded iCAP maps, in red, are displayed on the PAM50 template (De Leener et al., 2018) . Fine-grained iCAPs are overlaid over the corresponding atlas region probabilistic maps, in blue (Cadotte et al., 2015) . White ' + ' indicates the presence of bilateral clusters in one iCAP. All iCAPs are shown from rostral to caudal components and iCAP numbers are indicated in the bottom left corners. L, left; R, right; D, dorsal; V, ventral.
parate spatiotemporal properties. In particular, tSNR appeared larger using the OVS sequence, compared to the ZOOMit implementations. This reduction is to be expected, as tSNR is negatively linked to the size of the FOV (reduced by the square root of the FOV reduction factor) (Samson et al., 2016) (see Table 1 for sequence parameters). Likewise, the ZOOMit sequence with a slice thickness of 5 mm yielded a larger tSNR than for thinner 3 mm slices, as tSNR also decreases with voxel volume ( Murphy et al., 2007 ). On the other hand, sSNR followed an opposite trend, with both ZOOMit sequences exhibiting higher ratios than OVS, suggesting less image variation in the former case. We posit that this difference likely arises due to more severe signal losses when using OVS, at the level of intervertebral disks. It should be acknowledged that dedicated approaches have been proposed to address these issues, such as slice-specific shimming . However, they are For the three sequences, each fine-grained iCAP was matched to one of the 36 atlas regions defined by (Lévy et al., 2015) . We reported the number of iCAPs present in each region. Atlas regions 1 to 30 correspond to the white matter and regions 31 to 36 to the gray matter (separated by a dashed line). Odd (even) numbers are on the left (right) side. Dark (light) colors indicate regions found in iCAPs including unilateral (bilateral) cluster(s). Bilateral iCAPs are identified by different symbols in the middle of the corresponding bars (i.e., two clusters of the same iCAP identified by the same central symbol). not readily available in sequences provided by the manufacturer, thus hampering their deployment. SNR values were significantly higher in the gray matter than in the white matter, for all types of acquisition. While an opposite behavior, with larger values in the white matter, has previously been reported for cerebral tSNR measures ( Mazerolle et al., 2013 ), we believe that this discrepancy is likely related to the distinct organizational structure of the spinal cord, with the white matter being located outside the gray matter, at the interface with the CSF and closer to the venous plexus. Therefore, additional signal variations may arise in the white matter due to contribution from large vessels and CSF pulsations ( Piché et al., 2009 ). Importantly, our results are in line with earlier work performed in the spinal cord, where tSNR measures were also found to be higher in the gray matter .
Interestingly, characteristic motion patterns were observed for the different sequences. Specifically, motion was higher in runs performed with the OVS scheme. The consistency of this sequence dependent differences suggested their intrinsic nature. To better understand the origin of these fluctuations, we performed a spectral analysis of motion parameters. Our results revealed a peak in the range of the respiratory rhythm (ranging between 0.33 Hz and 1.66 Hz, (Fleming et al., 2011) , driven by translations in the y direction. A time-dependent voxel shift in the PE direction is, indeed, a typical signature of breathing-induced variations Fig. 7. iCAPs temporal properties. For the three sequences, we evaluated the temporal properties of the recovered fine-grained iCAPs. Mean values over iCAPs ± SE are presented for the total duration of summed activations ( A. ), the average duration of occurrences ( B. ), the total number of co-active iCAPs per time point ( C. ) and the couplings and anti-couplings, within-and between-level, as illustrated by Jaccard indices (i.e., joint activation time) ( D. ). In each plot, different sequences are represented by bars of different colors. * * corresponds to p < 0.01 and * * * to p < 0.001). ( Colloca et al., 2017 ;Verma and Cohen-Adad, 2014 ). This increased influence of respiration when using OVS might be associated with the larger FOV, leading to the inclusion of artifactual signal variations. Importantly, the combination of motion correction, motion scrubbing and physiological noise correction should effectively remove these artifacts ( Kassinopoulos and Mitsis, 2019 ).
In spite of these differences, all sequences demonstrated comparable performance in detecting motor-related activity during a wrist adduction movement. This advocates for the efficacy of our processing pipeline in reliably removing motion and physiological artifacts. Based on the innervation site of wrist adductor muscles ( Kendall et al., 2010 ), we hypothesized that the activity would be localized around the C7 and C8 spinal levels. In an earlier study combining muscle and spinal recordings (Kinany et al., 2019) , we have confirmed the distribution of the spinal pattern associated with this movement. While it is possible to probe the rostrocaudal profile of the activity obtained with the different acquisition schemes, it should however be mentioned that the type of movements used here (dynamic, bilateral and involving multiple muscles) elicits a combination of motor and sensory activity, thus hindering the precise in-plane comparison of the activation patterns. As expected, all three sequences exhibited a peak in the C7 spinal level, in line with our previous results. In a series of control analyses, we confirmed that spinal activation patterns were specific to task runs, further corroborating their motor origin. It should be emphasized that the statistical procedure used for GLM-based fMRI inference (correction for multiple comparison using a cluster-defining threshold of p < 0.001) is on par with brain standards. Consequently, the risk of inflated false positives is limited ( Eklund et al., 2016 ). Although group-level results obtained with OVS, Z3mm and Z5mm were similar, several observations suggested an increased sensitivity of OVS to task-related activity. This was notably reflected by the control analyses (larger activation patterns at the subject-level, along with higher Z scores), as well as by the superior repeat reliability. This may reflect the disparity in tSNR between the sequences, implying that a larger amount of data is needed to robustly image activation patterns using Z3mm and Z5mm.
The study of spontaneous activity is also of utmost interest, as it can clarify the spinal mechanisms supporting our body's physiological needs, as well as expose the functional backbone of spinal cord's intrinsic architecture ( Barry et al., 2014 ;Kinany et al., 2020 ;Kong et al., 2014 ). Additionally, resting-state recordings are ideal candidates for clinical applications, by means of their simplicity ( Fox, 2010 ;Lee et al., 2013 ). Considering this potential, we proceeded to compare the three acquisition schemes in this context. By deploying the SpiCiCAP framework ( Karahano ğlu and Van De Ville, 2015 ; Kinany et al., 2020 ), we showed that meaningful low-and high-granularity components, lining up with spinal levels and atlas regions, could be uncovered with all sequences. This supports our earlier results ( Kinany et al., 2020 ), advocating for the robustness and generalizability of the SpiCiCAP framework. Notwithstanding this, patterns obtained using ZOOMit sequences were more stable across subjects, hinting at a potential benefit of this approach to delineate spinal resting-state components. Another observation was the presence of bilateral iCAPs, mostly observed using OVS. Although we did not detect such patterns in our earlier study, using Z3mm ( Kinany et al., 2020 ), bilateral spinal resting-state networks have previously been reported in the literature ( Kong et al., 2014 ;Barry et al., 2014 ;Eippert et al., 2017b .). The prominence of these patterns in OVS might pertain to the lower sSNR for this sequence, possibly hampering its ability to disentangle neighboring regions. When evaluating the temporal properties of iCAPs derived from the three datasets, we differences in the reported metrics, as both the iCAP durations and the number of co-active iCAPs appeared to be dependent on the sequence. These differences may arise from the disparate TR values employed in the sequences. Nevertheless, a closer examination of temporal overlap underlined interaction patterns resembling the ones reported in our previous work ( Kinany et al., 2020 ). In particular, couplings were stronger than anti-couplings and most prominent within the same spinal level, alluding to the importance of local processing. Importantly, these patterns were also robust across acquisition schemes.
A conspicuous observation was the presence of alternative tracts, mostly using OVS and Z5mm, that were not detected in our earlier study ( Kinany et al., 2020 ). While components articulated around the CST and DCML pathways were observed, confirming our previous results, additional regions including ascending (spinocerebellar and spino-olivary) and descending (reticulospinal, rubrospinal, vestibulospinal and tectospinal) tracts were also exposed. These pathways differ from the CST and DCML pathways in several ways. First, they have smaller sizes and tend to be more scattered and intermingled ( Nathan et al., 1996 ). This makes them harder to distinguish and may have limited their detection. Then, they do not involve cortical regions, but are instead related to subcortical areas, such as the brainstem and the cerebellum (Rea, 2015) . Consequently, they carry different functions. While the CST and DCML are used for voluntary motor control and conscious sensory feedback, these alternative pathways mostly transmit unconscious information (e.g., unconscious proprioception) and support involuntary and automatic motor control (e.g., balance and posture). Although only conjectural, it is possible that the shorter TR of OVS and Z5mm (1.89 s and 1.5 s, respectively) have facilitated the detection of these regions. Further investigation is needed to better understand their origin and explore their spatiotemporal properties, though out of the scope of this comparative study.
Together, our findings indicate that the three acquisition schemes could be used to efficiently image spinal cord activity during task and at rest. An important contribution of this study was to highlight the robustness of spinal cord fMRI analyses against variations in the acquisition parameters, which advocates for the generalizability of spinal cord fMRI results. In spite of these overall similar performances, several observations could be drawn regarding sequence-specific properties, which translated into distinct advantages and drawbacks. Notably, the current study suggests that OVS may provide a superior sensitivity to task-related activation, while ZOOMit approaches may be more suitable to disentangle resting-state fluctuations. However, it is noteworthy that the choice of a particular set of parameters also depends on the researcher's needs and interests, as the experimental paradigm and research question will impact the technical requirements (e.g., coverage, spatial and temporal resolution, etc.). In line with recent efforts in establishing a reproducible pipeline for quantitative structural MRI (Cohen-Adad et al., 2021a, 2021b , next steps should include multicentric evaluations of spinal cord fMRI, so as to achieve standardized multi-vendor protocols for reliable assessments of spinal cord function. Although we used the ZOOMit sequence (Siemens), it should also be noted that similar implementations are provided by other manufacturers (e.g., FOCUS for GE and iZOOM for Philips), thus facilitating the deployment of this approach in different research and clinical centers. Finally, such spinal cord fMRI approaches could be translated beyond the cervical cord, to image other parts of the spinal cord, such as the thoracic or lumbar regions.
Importantly, the current study focused on investigating spinal cord fMRI signals acquired using specific sets of parameters, building on sequences previously deployed in the literature (e.g., Eippert et al., 2017b;Kinany et al., 2019;Weber et al., 2016 ). For this reason, sequences differ in multiple parameters (e.g., FOV, TR, TE, flip angle, etc.), which can, independently or interrelatedly, impact signal properties. As a result, while we demonstrated the robustness of the results for the three chosen acquisition schemes, we could not evaluate the role of each individual parameter. Future work is warranted to systematically probe the impact of selected parameters for a particular acquisition scheme (similarly to Barry et al. 2018 , for instance). Likewise, efforts are also required to assess to what extent task-and rest-related activity are sensitive to variations in the processing pipeline, expanding on earlier studies that have started undertaking this endeavor (e.g., Dehghani et al., 2020 for motion correction; Weber et al., 2017 for smoothing; Kong et al. 2012, Eippert et al., 2017a for denoising).
In conclusion, we compared three acquisition schemes for spinal cord fMRI and demonstrated their ability to reliably detect spinal cord activity during task and at rest. Our findings provide an auspicious indication of the generalizability of spinal cord fMRI results obtained using different acquisition parameters. Besides, we highlighted specific strengths and weaknesses associated with these different approaches to image spinal cord activity in vivo . Considering the unequivocal importance of the spinal cord in the human sensorimotor system, these methods can strongly benefit fundamental and clinical neuroscience.