FASB: an integrated processing pipeline for Functional Analysis of simultaneous Spinal cord-Brain fMRI

FASB: an integrated processing pipeline for Functional Analysis of simultaneous Spinal cord-Brain fMRI | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article FASB: an integrated processing pipeline for Functional Analysis of simultaneous Spinal cord-Brain fMRI Shahabeddin Vahdat, Caroline Landelle, Ovidiu Lungu, Benjamin De Leener, and 2 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3889284/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Simultaneous functional magnetic resonance imaging (fMRI) of the spinal cord and brain represents a powerful method for examining both ascending sensory and descending motor pathways in humans in vivo . However, its image acquisition protocols, and processing pipeline are less well established. This limitation is mainly due to technical difficulties related to spinal cord fMRI, and problems with the logistics stemming from a large field of view covering both brain and cervical cord. Here, we propose an acquisition protocol optimized for both anatomical and functional images, as well as an optimized integrated image processing pipeline, which consists of a novel approach for automatic modeling and mitigating the negative impact of spinal voxels with low temporal signal to noise ratio (tSNR). We validate our integrated pipeline, named FASB, using simultaneous fMRI data acquired during the performance of a motor task, as well as during resting-state conditions. We demonstrate that FASB outperforms the current spinal fMRI processing methods in three domains, including motion correction, registration to the spinal cord template, and improved detection power of the group-level analysis by removing the effects of participant-specific low tSNR voxels, typically observed at the disk level. Using FASB, we identify significant task-based activations in the expected sensorimotor network associated with a unilateral handgrip force production task across the entire central nervous system, including the contralateral sensorimotor cortex, thalamus, striatum, cerebellum, brainstem, as well as ipsilateral ventral horn at C5-C8 cervical levels. Additionally, our results show significant task-based functional connectivity between the key sensory and motor brain areas and the dorsal and ventral horns of the cervical cord. Overall, our proposed acquisition protocol and processing pipeline provide a robust method for characterizing the activation and functional connectivity of distinct cortical, subcortical, brainstem and spinal cord regions in humans. Biological sciences/Biological techniques/Imaging/Functional magnetic resonance imaging Biological sciences/Neuroscience/Neural circuit Biological sciences/Neuroscience/Sensorimotor processing Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Introduction A plethora of evidence suggests that many of the movement disorders, such as amyotrophic lateral sclerosis, multiple sclerosis, Parkinson’s disease, stroke, spinal cord injury, dystonia, and cerebral palsy impair both the brain and spinal cord circuits (Pierrot-Deseilligny and Burke, 2012 ; Raudino and Leva, 2012 ; Wahl and Schwab, 2014 ). However, human neuroimaging studies with patients affected by these disorders have mainly focused on characterizing the brain circuit dysfunction and deficits, while disregarding the critical role of the spinal cord circuitry and brain-spinal cord pathways. Simultaneous functional magnetic resonance imaging (fMRI) of the spinal cord and brain provides an unprecedented opportunity to address this limitation by allowing the assessment of the ascending sensory and descending motor pathways functioning in humans in vivo (Kinany et al., 2022 ; Landelle et al., 2021 ; Tinnermann et al., 2021 ). By employing this combined approach, researchers can investigate not only the neural activations at multiple levels of the central nervous system (CNS), including cortical, subcortical, brainstem, and spinal cord, but also functional interactions among these structures, thereby gaining valuable insights into the neural mechanisms underlying motor and sensory processing. Nevertheless, the imaging parameters and processing pipeline associated with this simultaneous neuroimaging approach are still in the early stages of development. In this study, we introduce an optimized acquisition protocol including the required anatomical and functional images, as well as an integrated analysis pipeline, named F unctional A nalysis of simultaneous S pinal cord- B rain fMRI (FASB). The proposed pipeline aims to improve the current challenges associated with the processing of spinal cord fMRI data under simultaneous spinal cord-brain acquisition in three domains. These improvements include motion correction, registration to the template, and removing the effects of large susceptibility artifacts, often exacerbated at the spinal disk levels. Prior work using simultaneous spinal cord-brain fMRI has revealed that task- and context-dependent dynamical interactions between brain regions and the spinal cord subserve various sensory and motor functions. In the last decade, the optimization of acquisition protocols at high-field strength (Barry et al., 2018b ), as well as the development of multi-slab pulse sequences designed to improve the acquisition of selective fields-of-view (FOVs), combined with advanced shimming procedures (Finsterbusch et al., 2013 , 2012 ; Islam et al., 2019 ) have significantly enhanced the quality of spinal cord-brain fMRI. Using these advanced image acquisition and shimming techniques, we have recently demonstrated the existence of integrated spinal cord and brain resting-state networks that follow neuroanatomical principles governing the ascending and descending pathways in humans (Vahdat et al., 2020 ). Despite clear advancement in spinal cord fMRI thanks to these new techniques, there is still work needed to improve image quality due to all the challenges associated with functional imaging of the spinal cord (Landelle et al., 2021 ). The application of spinal cord fMRI is starting to augment our brain-centered view of neurological disorders, by investigating the role of spinal cord circuit function in movement disorders, e.g. in Parkinson’s disease (Landelle et al., 2023 ). Yet studying the function of ascending and descending pathways requires the application of simultaneous fMRI methods, which has been limited in clinical populations due to a number of technical challenges (Landelle et al., 2021 ). These include the complex acquisition procedures that can achieve both a reasonable signal to noise ratio (SNR), and a high spatial resolution for distinguishing spinal cord quadrants, further constrained by the scanning time considerations. The goal of this study is to propose a simple acquisition method that can be implemented in clinical populations, together with an integrated processing pipeline to standardize the procedures for simultaneous spinal cord-brain fMRI. The spinal cord toolbox (SCT) (De Leener et al., 2017 ) has offered optimized set of tools for preprocessing of the spinal cord MRI data, including slice-wise motion correction, automatic segmentation of the spinal cord using convolutional neural networks (Gros et al., 2019 ), and nonlinear slice-wise registration to the template (De Leener et al., 2018 ). However, the combined acquisition of the spinal cord and brain provides new opportunities, as well as challenges, which can benefit from a customized processing pipeline. For example, the motion correction procedure implemented in the SCT is solely based on the spinal cord FOV, and can thus be improved by incorporating information from the acquired brain FOV. Also, large spatial distortions due to exacerbated susceptibility artifacts at the disk levels in BOLD imaging limit the efficacy of non-linear transformations for registration. Importantly, this may cause signal loss in the EPI data and may degrade SNR in some of the spinal cord slices. These low tSNR voxels practically degrade the group-level results by adding noisy voxels to the group-level averaging. Yet, there are two solutions to the latter problem. One approach is to remove the low SNR voxels from the participant’s spinal cord mask. However, the voxels included at the group-level analysis should be present in the mask of every single individuals. Because low-SNR voxels appear at different locations across participants (depending on neck curvature, surrounding tissues, etc.), a large number of voxels will be removed from the group-level analysis, which makes this approach impractical. Another approach is to remove the effect of these low-SNR voxels statistically from the group-level analysis, by incorporating participant-specific voxel-wise confounds in the group-level regression model. In this study, we employ the latter idea and propose a new approach that features automated modeling and identification of spinal cord voxels with low temporal SNR in each participant, followed by removing their effects statistically at the group-level analysis without the need to exclude those voxels from the entire group. To evaluate the efficacy of the proposed pipeline, FASB, we collected simultaneous spinal cord-brain fMRI data during the performance of a motor control task (unilateral handgrip force control), as well as during resting-state conditions. First, we compared the performance of the proposed algorithm for motion correction, registration to the spinal cord template, and physiological noise removal, to the current methods used in the spinal cord fMRI analysis. We then reported the activation maps in the brain, brainstem and spinal cord using FASB and other analysis methods during the motor control task. We hypothesized that the subject’s performance of this unilateral handgrip task would significantly increase the level of activity in the ipsilateral ventral horn at the C6-C8 cervical levels (Pierrot-Deseilligny and Burke, 2012 ), as well as in the contralateral primary and secondary sensorimotor cortical areas, striatum, thalamus, ipsilateral cerebellar cortex (Doyon et al., 2018 ; Vahdat et al., 2011 ), and ipsilateral motor-related brainstem nuclei, including pontine and medullary reticular formation (Arber and Costa, 2022 ; Ruder and Arber, 2019 ). We further assessed the functional connectivity of the cervical dorsal and ventral horns with various sensory and motor areas of the brain and brainstem. Our results demonstrate that the proposed acquisition protocol and processing pipeline (FASB) outperforms other methods commonly used in spinal cord fMRI analysis, and provides a robust method for accurately characterizing the activation and functional connectivity of distinct cortical, subcortical, brainstem and spinal cord regions in humans in vivo. Methods Participants This study was reviewed and approved by the ethics committee at the University of Florida, which adhered to the Declaration of Helsinki. All participants signed an informed consent form prior to participating in the study and were debriefed and compensated at the end of the experiment. Fifteen young, right-handed, healthy adults (10 females, mean age = 22.7 years, SD = 4.0) were selected to take part in this study based on the following exclusion criteria: presence of a history of neurological and psychiatric diseases, of any motor-system complication, the current use of any medication, and presence of any MRI-incompatible object in the body. An additional twenty-one healthy participants were recruited for resting-state acquisition (12 females, mean age = 31.4 years, SD = 7.1). These participants fulfilled the inclusion criteria previously described and the data acquisition was approved by the ethics committee at the University of Montréal, which also adheres to the Declaration of Helsinki. Due to poor fMRI image quality, resting-state data from two participants from this sample were not included in analyses, and data from another participant was also excluded due to missing physiological noise measurements. Thus, the resting-state sample included in the final analyses consisted of 18 healthy participants (mean age 32.39 (SD = 13.19), F/M = 11/7). Motor control task Participants who performed a motor control task in the scanner were asked to produce isometric hand grip force with their right hand using an MR-compatible fiber-optic force transducer. Prior to the MRI scan, participants were first trained on the task, and the grip maximum voluntary contraction (MVC) was measured. Specifically, they were asked to produce unilateral grip force using their right hand by opposing their thumb, index, and middle fingers against the force transducer (Archer et al., 2018 , 2016 ; Burciu et al., 2016 ). Two runs were completed in a block-design paradigm. Each run consisted of 9 blocks of 15 sec force control, followed by 25 sec rest periods (6 min and 25 sec in total). Throughout the run, 2 bars were displayed on a video monitor that participants saw through a mirror mounted on the head coil: a target bar and a force bar. A change in color of the force bar was used to cue participants to either push or release the force sensor. Green was a go signal for producing and sustaining force, while red indicated a rest period. The target bar was set at 10, 20, or 30% of MVC for each participant (3 blocks at each level, pseudo-randomly ordered), while the force bar moved in the vertical plane according to the force output. During the go signal, participants were required to adjust their grip force levels to match the force and target bars. Force was measured using custom-designed load cells (Vaillancourt et al., 2003 ) and a Micron Optics fiber optic interrogator. MRI data acquisition Motor control task dataset . A 3T Siemens Prisma scanner equipped with a 64-channel head/neck coil was used for imaging subjects. We acquired functional magnetic resonance imaging (fMRI) data, covering the whole brain and cervical spinal cord (first cervical, C1, to first thoracic, T1, vertebrae), using blood oxygenation level dependent (BOLD) contrast via a simultaneous multi-slice (SMS) echo-planar imaging (EPI) sequence. Participants laid on the scanner table in a supine position with their head and neck fully supported using foam pads to minimize their motion. We used the neck and brachial plexus SatPad™ sets (SatPad, clinical imaging solutions) around the neck to minimize inhomogeneity artifacts and enhance spinal cord functional image quality. Based on previous work in the lab, we recommend their use in fMRI studies that target the spinal cord. Participants were placed in the scanner such that the mid-chin level was in the scanner’s isocenter, corresponding to C2-C3 vertebral level on the axial plane. The following parameters were used for EPI acquisitions: 72 axial slices angled orthogonal to the cord at C5 vertebral level, in-plane resolution: 1.6x1.6 mm 2 , matrix size: 120x120, 50% phase oversampling, slice thickness: 4 mm (no gap), TR: 1.85 s, 210 volumes, TE: 27 ms, multiband acceleration factor: 3, GRAPPA factor: 2, 0.875 partial Fourier encoding, FA: 70˚, PE direction: A-P. Additionally, 2 saturation pulses were applied in a V-shaped configuration to minimize ghosting and inflow artifacts related to blood flow in the major cervical vessels. A 3D-MPRAGE sequence was used to acquire T1-weighted anatomical images covering the entire head as well as neck down to the T3 vertebral level using the following parameters: field of view (FOV): 192×240×320 mm 3 , sagittal slices; TR: 2300 ms, TE: 3.41 ms, and resolution: 1×1×1 mm 3 . In addition, a multi-echo data image combination (MEDIC) sequence was used to acquire a T2*-weighted (T2w) anatomical image using the following parameters: axial slices = 36, slice thickness = 4 mm (no gap), in-plane resolution: 0.75×0.75 mm 2 , matrix size: 256×256, TR: 1280 ms, TE: 18, 20, 21, and 22 ms; FA: 20˚, and GRAPPA acceleration factor: 2. Although the orientation and the size of slices were matched between the T2w and EPI images, the T2w image had fewer slices, covering only the cervical cord to keep the acquisition time down to 12 min. Importantly, the position of the T2w image FOV in the dorsocaudal orientation was set such that the position of the most caudal slice matches that of the EPI image. Resting-state dataset. Data were also acquired using a 3-Tesla MRI Scanner (Magnetom-Prisma, Siemens, Erlangen, Germany) equipped with a 64-channel head (1–7 elements active) and neck coil (1–2 elements active). T2-weighted images were collected in the transverse orientation with the following parameters: TR = 33ms; TE = 14ms; FOV = 224 mm x 258 mm; flip angle = 5 o ; in-plane voxel resolution = 0.35 x 0.35 mm 2 , slice thickness = 2 mm. A total of 112 slices were acquired, spanning from the top of the cerebellum to the upper thoracic region (approximately at T1), thereby encompassing the entirety of the cervical spinal cord. Resting-state functional images were acquired using an echo-planar imaging (EPI) gradient echo sequence covering the brain and cervical spinal (TR/TE = 1550/23 ms, axial FOV = 120 x 120 mm 2 , flip angle = 70°, in-plane voxel resolution = 1.6 x 1.6 mm 2 ; slice thickness = 4 mm (no gap); iPAT acceleration factor PE = 2, iPAT acceleration factor slice = 3). A total of 69 axial slices were acquired covering the top of the head to approximately T1 vertebra. Only the slices containing the cervical spinal cord were included in further analyses. One resting state-run ( i.e. , no explicit task) was acquired in total with an acquisition time of 10 minutes and 35 seconds each (230 volumes). Participants were instructed to relax watch the video and minimize swallowing and motion. Physiological recordings were acquired using a pulse sensor and a respiration belt (Siemens Physiology Monitoring Unit). Spinal cord analysis pipeline We utilized bash programming, the FMRIB Software Library (FSL) (Smith et al., 2004 ), Spinal Cord Toolbox (SCT version 6.0) (De Leener et al., 2017 ), and in-house MATLAB programs to process the collected data. The full spinal cord processing pipeline will be publicly available at www.vahdatlab.org . The main steps of the spinal cord processing pipeline are shown in Fig. 1 , including preprocessing, registration to the spinal cord template (PAM50 from SCT toolbox), and regression analysis, as described below. Preprocessing Following removing of the first few volumes to reach a steady state condition (in our dataset, no volume was removed due to the inclusion of dummy scans during the acquisition), motion correction was performed on the EPI data in two steps (Fig. 1 ). First, a 3d linear rigid-body transformation was performed using mcflirt tool (FSL) on the entire brain and spinal cord FOV. Next, the EPI was divided into the spinal cord and brain FOVs. Slice-wise (2d) non-linear motion correction was then performed in a relatively large (23 voxels diameter) cylindrical mask of the spinal cord using the sct_create_mask (SCT). The performance of various motion correction methods (3d linear alone, slice-wise nonlinear alone, or both) was calculated and compared using two indices: i) DVARS that measures changes in BOLD signal intensity from one volume to the next (Eq. 1, (Afyouni and Nichols, 2018 ), and ii) tSNR (temporal signal-to-noise-ratio) which is the ratio of the mean by the standard deviation of the BOLD signal timeseries ( \({Y}_{t}\) ) in each voxel (Eq. 2), n being the number of acquired volumes. Equation 1 \(\text{D}\text{V}\text{A}\text{R}\text{S}=\sqrt{\frac{1}{n}\sum _{t=1}^{n}{({Y}_{t}-{Y}_{t-1})}^{2}}\) Equation 2 \(\text{t}\text{S}\text{N}\text{R}= mean\left({Y}_{t}\right)/std\left({Y}_{t}\right)\) Following motion correction, gradient-echo field map correction was performed using FUGUE tool (FSL) to correct for geometrical distortion due to magnetic field inhomogeneities (Fig. 1 ). These artifacts are usually larger in the spinal cord compared to the brain due to the proximity to the thorax, vertebrae, and disks, and can cause a shift in the phase-encoding direction, which can be compensated at this step. Next, spatial smoothing was performed inside a mask of the spinal cord in a slice-wise fashion using 3dBlurToFWHM command from AFNI (Cox, 1996 ). A kernel size of 3.2 mm (2 times the voxel size) was used for smoothing, although a 2.4 mm kernel size resulted in very similar task-based activation maps. Registration to template The spinal cord registration pipeline is displayed in Fig. 2 . Prior to registration, a spinal cord mask was generated for each of the T2w and T1w images using a deep learning segmentation algorithm as implemented in sct_deepseg_sc (SCT) with automatic centerline detection using the support vector machine and convolutional neural network options (Gros et al., 2019 ). Following visual inspection, manual corrections were carried out to address any minor inaccuracies in the spinal cord mask. Next, a centerline was drawn manually on the T2w and EPI images. To be consistent across participants, in the T2w image, the centerline in the dorsoventral direction was defined at a voxel location just posterior to the gray commissure, which was clearly visible in the acquired T2w image. For the EPI image, since the gray/white matter distinction is less clear, we used the calculated tSNR map to define the centerline. As shown in Fig. 3 , the cord exhibits noticeably higher tSNR values compared to the surroundings, including the cerebrospinal fluid (CSF), blood vessels, and connective tissues. Hence, we identified the centerline at the midpoint of the area with the highest tSNR (yellow area in Fig. 3 ; when tSNR map is thresholded at 40, color-coded in red-yellow). Registration to the template (PAM50) was performed in three steps (Fig. 2 ): 1) EPI image registration to the T2w space using slice-wise centerline images alignment, 2) T2w image registration to the T1w space through the cord segmentation images using sct_register_multimodal (SCT), and 3) T1w image registration to PAM50 through cord segmentation using sct_register_to_template (SCT). In detail, for step 1, EPI and T2w images were registered by performing slice-wise 2d in-plane translations (in the dorsoventral and mediolateral directions) to align their centerlines using flirt command (FSL) with the nearest-neighbor interpolation. For step 2, we manually identified the C1-T1 disc labels in both T1w and T2w images, and used their calculated cord segmentations images to compute a slice-wise rigid body transformation using sct_register_multimodal command (any mismatch in the rostrocaudal direction was first corrected by aligning T1w and T2w disk label locations). For step 3, the cord segmentation and disc labels in the T1w image were employed for registration to PAM50 through sct_register_to_template command, using the default options (including straightening, vertebral alignment, and non-linear slice-wise registration). Note that only step 1 was applied to the 4d EPI images. The steps 2 and 3 warping transformations were concatenated using sct_concat_transfo , and were later applied to the summary statistics images (regression coefficients and their variance images) using sct_apply_transfo command . Subsequently, we compared the performance of the proposed registration method with the gold-standard spinal cord registration approach which only utilizes the T1W image (Landelle et al., 2023 ; Vahdat et al., 2015 ). This approach was implemented by calculating the transformations from the EPI to T1W image using sct_register_multimodal , and from T1W image to PAM50 as calculated above. The two transformations are then concatenated to register the EPI image directly to PAM50. To evaluate the registration performance, the binary cord segmentation mask in the EPI space was registered to PAM50 space using each method, and the result was compared to the binary template cord mask. Two indices were defined by calculating: i) the number of voxels in the template mask that are missed in the registered mask, normalized by the total number of template cord voxels, and ii) the overlap between the registered and the template cord masks, normalized by the total number of voxels in both masks. Regression analysis We performed participant-level general linear model (GLM) analyses on the preprocessed EPI data in the T2w space in two steps. First, a GLM was constructed to model and remove the effects of three sets of confounds regressors, including: i) the mean CSF signal averaged in a surround mask generated by dilating the mask of spinal cord ( i.e. , GM + WM), ii) 6 motion correction parameters (x, y, and z translations and rotations) derived from the linear motion correction step, and 2 slice-wise motion correction regressors generated by the non-linear motion correction step on the spinal cord FOV, and iii) 5 regressors, identified as confound, extracted from independent component analysis (ICA). The fastica algorithm (Hyvarinen, 1999 ) was used to decompose the preprocessed spinal cord EPI data (before spatial smoothing) within a large mask (including both the spinal cord and the surround mask) into 30 components (Vahdat et al., 2020 ). We calculated 3 indicators to systematically sort components based on their likelihood of being task-related, as opposed to those related to physiological/scanner noise: 1) spatial indicator: the ratio of significantly active ( p < 0.05) voxels within the spinal cord mask to that within the surround mask, 2) frequency indicator: the ratio of the power of the component’s time-series frequency response around the task central frequency (in our block-design paradigm: 1 / (task period + rest period) = 1/40 s = 0.025 Hz; 0.01 < f < 0.1 Hz may be used for resting-state data) to that outside the task-related frequency range with some margins ( f 0.025*1.5 Hz), and 3) temporal indicator: the Pearson correlation between the task timing signal (boxcar function in our task) convolved with a canonical hemodynamic response function (HRF) and the time-series of each component. A confound index was defined by multiplying the 3 indicators; the lower the confound index , the higher the likelihood that component represents a confound. The 5 components with the least confound index were selected, and their time series were included in the first GLM as confound. Next the residual (clean) image from the first GLM was high pass filtered at 1/60 Hz. The timing of the grip force block design task was modeled as a regressor of interest, while motion outliers were modeled as confounds. The latter were calculated using fsl_motion_outliers command on the non-preprocessed spinal cord EPI data (Fig. 1 ). The summary statistics images ( Cope and VarCope images in FSL) related to the regressor of interest were then registered to PAM50 space using the concatenated warping transformation calculated in the registration step. We generated a set of individual-specific masks to specify low tSNR voxels for each participant as follows. For each EPI run, a tSNR threshold was calculated using Maximum a Posteriori (MAP) estimation of two Gaussian distributions on the calculated tSNR histogram data inside a large spinal mask (including both the cord and the surround binary masks; Fig. 3 ). A Gaussian mixture distribution model using two components was fitted to the tSNR histogram using MATLAB fitgmdist function. The threshold was defined as the MAP estimate of the optimal boundary, which is the intersection of the two Gaussian distributions (Fig. 3 ). Low tSNR voxels (tSNR < threshold ) were identified for every EPI run and participant. A subject-specific low tSNR mask was then defined as the intersection of low tSNR voxels from all EPI runs for each participant (a total of 2 runs in our data). Note that the two Gaussian distributions represent, hypothetically, the spinal cord and the CSF/surround voxels, with high tSNR voxels expected to be located inside the spinal cord mask as shown in Fig. 3 . Hence, the goal of this step is to calculate an optimal threshold for differentiating spinal cord and CSF/surround voxels and identify those voxels inside the spinal cord mask that, inconsistently, possess low tSNR due to severe distortion and signal loss in the EPI spinal data. Finally, a group-level GLM model was constructed to average individual-level maps using a mixed-effect (FLAME 1 + 2) model as implemented in FSL. The individual-specific low tSNR masks modeling the low tSNR voxels in each participant were entered as confound (Fig. 1 ). First we generated a group-level spinal cord mask that contains voxels, in which the majority of participants showed an acceptable tSNR ( i.e. , tSNR > threshold ). This condition ensures that for every voxel in the group-level mask, high quality data from at least half of the participants were included for averaging, and, at the same time, limits the number of confound regressors to half the degrees of freedom in the group-level GLM. Hence, a total of up to n = (number of participants / 2) 4d voxel-wise regressors were included to account for and remove the effect of low tSNR voxels in every participants from group-level averaging. This procedure enabled us to only incorporate voxels with an acceptable tSNR in the group-level averaging, without the need to remove these voxels from the group-level mask/analysis. All group-level statistical maps were then corrected for multiple comparisons using Gaussian random field (GRF) theory, Z > 2.7, and cluster-level threshold p < 0.05. Brain fMRI processing Brain image preprocessing was performed with the FSL software package, using the same pipeline as described previously (Vahdat et al., 2020 , 2015 ). One participant was removed from brain analysis due to partial coverage of brain FOV. Preprocessing included the following steps: skull stripping (FSL, BET), motion correction (FSL, MCFLIRT), high-pass temporal filtering (1/60 Hz), GRE field map correction, spatial smoothing using an isotropic Gaussian kernel of 3.2 mm full width at half maximum. BBR method (FSL) was used to register the EPI data to the T1w anatomical image, followed by a nonlinear registration (FNIRT, FSL) from the T1w image to the MNI template. These two transformations were concatenated to register the summary statistics images from the participant level analysis in the EPI space to the MNI template. Participant-level GLM was constructed to identify brain areas activated during motor control task using a boxcar function modeling our block-design paradigm. The 6 motion correction parameters (x, y, and z translations and rotations) derived from the linear motion correction step were included as confound in this analysis. Furthermore, volumes detected as motion outliers (described earlier) were included as confound in the GLM. Next, a group-level GLM model was constructed to average participant-level maps using a mixed-effect (FLAME 1 + 2) model as implemented in FSL. All group-level statistical maps were then corrected for multiple comparisons using Gaussian random field (GRF) theory, Z > 2.7 (voxel-wise uncorrected threshold), and a cluster-wise corrected threshold at p < 0.05. Functional connectivity analysis We performed region of interest (ROI)-based functional connectivity analysis to examine the communication between the spinal cord and brain BOLD time-series. We defined 11 a-priori selected ROIs in the brain shown to be part of the sensorimotor network using the MNI brain atlas and brainstem atlas (Cauzzo et al., 2022 ; Singh et al., 2021 ; Vahdat et al., 2020 ). The brain ROIs included: left precentral gyrus, left postcentral gyrus, supplementary motor area, left thalamus, right cerebellum lobule I-IV, right cerebellum lobule V-VI, left red nucleus, bilateral pontine reticular formation, bilateral medullary reticular formation. We also defined 2 ROIs in the spinal cord (right ventral horn C5-C8, right dorsal horn C5-C8 segmental levels). To restrict our analysis to subregions somatotopically related to the forelimb region, each ROI only included the most active voxels during handgrip force production by selecting voxels with the highest Z -score from the individual-level GLM analysis (top 10, 20, 30, or 40th percentile voxels within each ROI). We then calculated the average clean BOLD signal (the residual image from the first GLM) within each ROI at each percentile threshold. For each participant and run, the partial correlation (Marrelec et al., 2006 ; Vahdat et al., 2014 ) was calculated between each pair of the brain and spinal cord ROI’s time-series (22 pairs), by accounting for and removing the effects of task presentation (task timing convolved with HRF). The average partial correlation across 4 percentile thresholds were calculated and reported as the functional connectivity (FC) for each ROI pair. For each participant, FC values between the two runs were averaged and Fisher transformation was applied to meet normal distribution assumption. Finally, t -statistics were performed to identify significant functional connections during the motor task between the brain and spial cord ROIs across participants (corrected for multiple comparisons using Bonferroni correction, p < 0.05). Resting-state analysis Denoising pipelines. To investigate the impact of different sources of noise we compared the performance of different denoising methods. For each methods, the fMRI time-series was regressed with various variables of no-interest (confounds) using FSL's fMRI Expert Analysis Tool (FEAT). No filtering was applied at this step. All methods incorporated 8 motion-related regressors, 6 obtained from the linear motion correction step, and 2 slice-wise regressors from the SCT motion correction step, as described earlier. Next, we extracted the mean signal of the CSF in a surround mask generated by dilating the mask of the spinal cord. For each participant, the physiological noise modelling was calculated using the Tapas Physio toolbox (Kasper et al., 2017 ). We used the RETROspective Image CORrection (RETROICOR) procedure (Glover et al., 2000 ) to generate noise regressors from peripheral physiological recordings (heart rate and respiration). Specifically, we modeled, four cardiac, four respiratory, harmonics, and four multiplicative terms for the interactions between respiratory and cardiac noise (32 regressors in total, similar to (Brooks et al., 2008 ; Kaptan et al., 2023 ; Kong et al., 2012 ). To identify non-neural fluctuations, we extracted the first five principal components of the unfiltered and unsmoothed CSF signal (in the surround mask) in the participant’s native space using the CompCor approach (Behzadi et al., 2007 ). Additional regressors, identified as confound, were extracted from ICA in a dilated mask of the spinal cord as described in the Regression analysis section (3.89 ± 1.66 components). Finally, a specific set of regressors of no-interest was accounted for in separate GLMs, corresponding to different denoising models as follow: “ MeanCSF ” (baseline parameters, mean CSF) “ ICA ” (baseline parameters, mean CSF, ICA components) “ CompCor ” (baseline parameters, mean CSF, 5 CompCor components) “ Retroicor ” (baseline parameters, mean CSF, 32 Retroicor components) “ ICA + CompCor ” (baseline parameters, mean CSF, ICA and 5 CompCor components) Variance of the model. For each method, we calculated the voxel-wise explained variance of the model (R 2 ) as a performance metric to provide deeper insights into the impact of noise removal from various sources using the following equation: Equation 3 \({\text{R}}^{2}=[var\left({Y}_{t}\right) - var\left({\epsilon }_{t}\right)]/var\left({Y}_{t}\right)\) Where \(var\left({Y}_{t}\right)\) and \(var\left({\epsilon }_{t}\right)\) denote respectively the variances of the fMRI times-series and the residual signal ( res4d.nii in FSL) at each voxel, after applying the regression model associated with each method. Next, the average R 2 was calculated inside the surround mask ( i.e. , the CSF mask), as well as the spinal cord mask. The ratio between the R 2 of CSF and spinal cord (R 2 CSF / R 2 cord ) is reported to evaluate the model's performance. A good model for removing the physiological noise would account for more variance in the CSF, while for less variance related to neural activation in the spinal cord. Hence, the larger the R 2 CSF / R 2 cord ratio, the greater the capacity of the model in explaining the BOLD signal variability inside the CSF compared to the spinal cord, and the better the model in removing the physiological noise. Being aware of the limitations of the R 2 as a quality check metric, such as its failure to account for each model's degree of freedom, we also calculated the adjusted variance of the model (AdjR 2 ) as follows: Equation 4 \({\text{A}\text{d}\text{j}\text{R}}^{2}=1- \frac{(1-{\text{R}}^{2})(N-1)}{N-P-1}\) Where N is the number of fMRI image time points, and P corresponds to the number of regressors of no-interests (confounds) included in each model. Unlike the R 2 , which measures the proportion of variance explained by the model, AdjR 2 takes into account the degree of freedom of the model. This implies that with an increase in the number of regressors (resulting in the degree of freedom increase), the AdjR 2 value will increase only if the additional regressors enhance the overall model performance. The AdjR 2 was extracted in the CSF mask. t-statistics were performed to compare the R 2 and AjdR 2 between each model, and differences were considered significant with a p -value < 0.005, which was adjusted for multiple comparisons using Bonferroni correction (0.05 /10 = 0.005). Resting-state functional connectivity. To assess the impact of denoising methods on functional connectivity measures, we conducted ROI-based analyses on the residual (cleaned) images. We defined 4 ROIs in the spinal cord (right ventral horn C5-C8, right dorsal horn C5-C8, left ventral horn C5-C8, left dorsal horn C5-C8) and 1 ROI in the surround mask cord (CSF C5-C8). We extracted and band-pass filtered (0.01–0.17 Hz)(Barry et al., 2018a ) the time series within each ROI and calculated the average signal. Next, for each participant, we computed the partial correlation between each pair of spinal cord ROIs while the CSF average signal was used as a covariate. The Fisher-transformed correlation coefficients were employed to establish functional connectivity between each ROI pair, specifically between ventral horns (Ventro-Ventral FC), dorsal horns (Dorso-Dorsal DC), ventro-dorsal right horns (Right FC) and ventro-dorsal left horns (Left FC). For each FC, t-statistics were performed to compare the functional connectivity between each model (corrected for multiple comparisons using Bonferroni correction 0.05 /10 = 0.005). Results We collected simultaneous spinal cord-brain fMRI data, covering the whole brain and the entire cervical cord down to the first thoracic vertebra, using simultaneous multi-slice (SMS) EPI sequence, together with the application of neck and brachial plexus SatPad™ set, as described in the Methods section. This image acquisition protocol provided a reasonable spatial and temporal resolution (1.6x1.6 mm 2 in plane, TR: 1.85 ms), as well as a reasonable tSNR in the spinal cord (22.32 ± 0.86 mean ± sem; Fig. 3 ) of our group of participants. Furthermore, the low tSNR voxels were statistically excluded, per participant, from the group-level analysis. On average 76 ± 2.1 (mean ± sem) percent of the spinal cord voxels passed the MAP threshold and were included in group-level analysis. The average tSNR in the voxels included in the group-level analysis was 25.40 ± 0.81 (mean ± sem) across participants (Fig. 3 ). Motion correction procedures First, we assessed the efficiency of various motion correction procedures in the spinal cord using DVARS (Eq. 1) and the spinal cord tSNR metrics (Eq. 2). Optimal outcomes encompass reduced DVARS and enhanced spinal cord tSNR values. The performance of 3 methods were compared to the raw (no motion correction) data, including 3d linear motion correction on the entire brain and spinal cord FOV (SP + BR Flirt), slice-wise nonlinear correction on the spinal cord FOV (Slice-wise SCT), and the combined procedure (first performing 3d linear correction on the entire FOV, followed by slice-wise nonlinear correction on the spinal cord slices) as implemented in FASB. As shown in Fig. 4 , the combined procedure outperforms the other methods in terms of DVARS and tSNR (repeated-measures ANOVA, corrected for multiple comparisons using Bonferroni correction, DVARS : F(3,42) = 14.3, p (FASB vs Raw) < 0.005, p (FASB vs Flirt) < 0.05, p (FASB vs SCT) 0.05; Spinal cord tSNR : F(3,42) = 23.5, p (FASB vs Raw) < 0.001, p (FASB vs Flirt) < 0.01, p (FASB vs SCT) 0.05). This analysis revealed that combining both 3d linear and slice-wise non-linear motion correction procedures yields reduced temporal variability and enhanced temporal SNR in the spinal cord. Registration procedures Next, we compared the performance of various registration methods to align the EPI spinal cord image to a template (PAM50). As outcome measures, we calculated percentage error and overlap between the registered binary EPI spinal cord mask and the PAM50 spinal cord segmentation, as reported in Methods. Smaller Error and larger Overlap percentage values are desirable. As described in the Methods, the performances of 3 procedures are compared, including registration using T1W image as the intermediate image (T1W-based), GRE field map correction together with T1W-based registration (GRE + T1W-based), and FASB registration procedure as summarized in Fig. 2 . As shown in Fig. 5 , the FASB registration procedure outperforms the other registration methods (which do not include T2w image) in terms of both Error and Overlap for the spinal cord masks (repeated-measures ANOVA, corrected for multiple comparisons using Bonferroni correction, Error : F(2,28) = 37, p (FASB vs T1W) < 0.001, p (FASB vs GRE + T1W) 0.05; Overlap : F(2,28) = 11.5, p (FASB vs T1W) < 0.001, p (FASB vs GRE + T1W) 0.05). This analysis shows that performing the additional slice-wise centerline alignment between EPI and T2w images, as proposed in FASB, significantly improves the overall accuracy of the spinal cord registration to the template space (from 78–84% overlap). Physiological noise reduction procedures Next, we compared the performance of different physiological noise removal models on the collected resting-state fMRI data. We compared 5 models (see Methods for details): “meanCSF” (the baseline model for removing motion correction parameters and the CSF average), “ICA” (as implemented in FASB), “CompCor” (a PCA-based model), “Retroico r ” (based on the recorded cardiac and respiratory signals), and “ICA + CompCor” (combined “ICA” and “CompCor” models). First, we compared the ratio between the R 2 of CSF and spinal cord (R 2 CSF / R 2 cord ) between all denoising models to evaluate model performances (Fig. 6 A). A higher value indicates that the variance of the CSF is explained more substantially by the model relative the variance of the cord. Results showed a significant increase across most of the denoising models compared to the baseline model “ meanCSF ” (" ICA ” vs . “ meanCSF : t(17) = 4.67, p = 0.0002; " CompCor ” vs . “ meanCSF” : t(17) = 4.87, p = 0.0001; " ICA + CompCor ” vs . “ meanCSF : t(17) = 4.34, p = 0.0004), with the exception of the “ Retroico r” model where no significant difference was observed with the “ meanCSF ” model (t(17) = 1.12, p < 0.27 ). In addition, “ Retroico r” model showed a lower ratio than “CompCor” (t(17)= -5.38, p < 0.0001), “ICA” (t(17)= -4.02, p = 0.0009) and “ ICA + CompCor ” (t(17) = 5.97, p < 0.0001) models. “CompCor” model had significant higher ratio than “ICA” ( " CompCor ” vs . “ ICA” (t(17) = 3.45, p = 0.003) , while no significant differences, after the Bonferroni correction, were observed between the “ICA ” and “ ICA + CompCor” (t(17) = 3.03, p = 0.007 ) , or " CompCor ” and “ ICA + CompCor” (t(17) = 2.60, p = 0.018 ). To account for the degrees of freedom of each model (depending on the number of confound regressors), we also calculated the adjusted R 2 (AdjR 2 ; Eq. 4) inside the CSF mask. Adding new confound regressors does increase the AdjR 2 value only if the additional regressors enhance the overall model performance. As shown in Fig. 6 B, all models have higher AdjR 2 than the baseline model ( “ICA ” vs . “ meanCSF ”: t(17) = 12.89, p < 0.0001; “CompCor ” vs . “ meanCSF ” : t(17) = 24.25, p < 0.0001; “Retroicor ” vs . “ meanCSF”: t(17) = 7.01, p < 0.0001; “ICA + CompCor ” vs . “ meanCSF ”: t(17) = 24.38, p < 0.0001). “ Retroico r” model showed lower AdjR 2 than “ CompCor ” (t(17)= -20.40, p = 0.0001), “ICA” (t(17)= -7.07, p < 0.0001) and “ ICA + CompCor” (t(17)= -19.45, p < 0.0001) models. “ICA + CompCor” model had significantly higher AdjR 2 than “ICA” (t(17) = 20.40, p < 0.0001) and “CompCor” models (t(17) = 9.33, p < 0001). Finally “ CompCor” model showed significantly higher AdjR 2 than “ICA” model (t(17) = 10.35, p < 0001). Assessing the impact of denoising procedures on functional connectivity Next to assess the impact of denoising methods on functional connectivity (FC) measures, we conducted ROI-based analyses (see Methods). As shown in Fig. 6 C, no significant differences in FC between the different denoising models were observed associated with ipsilateral functional connectivity measures (ventro-dorsal right, and ventro-dorsal left). “ICA + CompCor” model led to significantly reduced functional connectivity between the left and right ventral horns compared to most of the other models (“ meanCSF ”: t(17)= -3.86, p = 0.001; “ ICA ”: t(17)= -3.80; p = 0.001; “ Retroicor ”: t(17)= -4.36, p = 0004) and a trend toward significant when compared to the “ CompCor ” model (“ ICA + CompCor ” vs . “ CompCor ” t(17) = 3.1, p = 0.006). Similarly, the FC between the left and right dorsal horns was significantly lower for the “ ICA + CompCor ” model compared to the others (“ meanCSF ”: t(17)= -5.59, p < 0.0001, “ ICA ”: t(17)= -4.02, p = 0.0008, “ Compcor ”: t(17)= -3.62, p = 0.002, “ Retroicor ”: t(17)= -6.19, p < 0.0001). The FC was also found lower with “ICA” model compared to “ meanCSF ” (t(17):-4.46, p = 0.0003) and “ Retroicor ” (t(17)= -4.56, p = 0.0002) models. In sum, ventro-ventral and dorso-dorsal connections showed reduction in functional connectivity in the more complex model (“ ICA + CompCor ”) compared to the others. Interestingly for the weaker FC observed within horns (ventro-dorsal FC), the various denoising models did not significantly impact the FC. Task-related results Figure 7 shows activation clusters associated with the motor task performance in the brain and brainstem. At the brain level, all of the anticipated sensorimotor cortical and subcortical areas related to unilateral hand movement (Vahdat et al., 2017 , 2015 , 2011 ) showed significant activation, including the contralateral primary motor cortex (M1), dorsal premotor cortex, primary and secondary somatosensory cortices (SI, SII), thalamus (ventral lateral and pulvinar nuclei), putamen, ipsilateral cerebellar cortex (lobules I-IV, and V-VI), bilateral posterior parietal cortex (PPC), and supplementary motor area (SMA) (corrected for multiple comparisons using GRF, Z > 2.7, p < 0.5, Fig. 7 A). Interestingly, our high in-plane spatial resolution (1.6×1.6 mm 2 ) allowed us to detect and delineate task-related brainstem activations, including bilateral pontine reticular formation (PRF), medullary reticular formation (MRF), and red nucleus (RN) (Fig. 7 B). Consistent with the laterality of descending projections from the reticular formation and RN (Davidson et al., 2007 ; Lemon, 2008 ; Soteropoulos et al., 2012 ), task-related activation was greater in the ipsilateral MRF and PRF, and the contralateral RN. These results show that the acquired EPI brain images in our simultaneous protocol not only provide enough power to replicate previous task-related brain networks, but also offer detection and delineation of activity in important brainstem nuclei. Next, we also compared task-related activation maps in the spinal cord using different approaches. Four methods were implemented and compared as described in the Methods section: i) regular group-level regression, without the inclusion of ICA confounds (NoICA-NotSNR), ii) regular group-level regression, with the inclusion of ICA confounds (ICA-NotSNR), iii) group-level regression by inclusion of voxel-wise low-tSNR confounds, but without ICA confounds (NoICA-tSNR), and iv) group-level regression by inclusion of voxel-wise low-tSNR confounds, as well as ICA confounds (ICA-tSNR, implemented in FASB). As shown in Fig. 8 , all 4 methods resulted in significant activation clusters at the C5-C8 cervical segmental levels (corrected for multiple comparisons using GRF, Z > 2.7, p < 0.05). Compared to the other methods, FASB (ICA-tSNR) resulted in the largest number of activation clusters (4 clusters, including C5, C6, C7, and C8) and detected activation across all related cervical levels for the control hand/wrist muscles (Landelle et al., 2021 ). Importantly, the largest peaks of activation were unequivocally located in the ipsilateral ventral horn, where the hand/wrist motoneurons are located. The activation at the C6-C7 level was largely restricted to the ventral horn, while it included the dorsal horn at the C8 segmental level. These results show that the acquisition/analysis pipeline proposed in FASB can differentiate task-related spinal cord activation based on myotome/dermatome (spinal level), side (left/right), and sensory-motor (dorsal/ventral) characteristics. To evaluate the connectivity related to the ascending and descending projections, we next investigated the partial correlation between a-priori selected set of ROIs in the brain and spinal cord, when the effects of task paradigm are removed from this relationship (see Methods for details). The right (ipsilateral) ventral horn at C5-C8 cervical level showed significant functional connectivity to the left postcentral gyrus, SMA, thalamus, and right medullary reticular formation and cerebellum lobule V-VI (Fig. 9 ; p < 0.05). Also, the right dorsal horn at C5-C8 cervical level showed significant functional connectivity to the left precentral and postcentral gyrus, SMA, thalamus, and right cerebellum lobule V-VI. These results demonstrate that significant task-related functional connectivity between the brain and spinal cord ROIs can be identified using FASB, consistent with the neuroanatomy of well-established descending and ascending projections in humans (Landelle et al., 2021 ; Lemon, 2008 ). Discussion We developed and validated an integrated processing pipeline for functional analysis of simultaneous spinal cord-brain EPI data (FASB), together with recommendations on collecting appropriate anatomical and functional data for that pipeline. The proposed EPI acquisition allowed us to acquire spinal cord images using a Siemens 3T MRI scanner, with low distortion due to magnetic susceptibility artifacts, and with high SNR (mean temporal SNR > 20 in the cord), while having a relatively high spatial resolution at the spinal cord level (1.6 mm in-plane), and a reasonable acquisition time (TR = 1.85 s). The novel aspects of the FASB pipeline are threefold (Fig. 1 ): 1) Improved motion correction of spinal cord fMRI volumes using both 3d linear and slice-wise non-linear transformations, 2) Improved spinal cord EPI registration to the template space (PAM50) using two intermediate anatomical images (T2w and T1w), and slice-wise centerline alignment to T2w image guided by the tSNR map of the EPI data (Fig. 2 ), and 3) Improved detection power of the group-level analysis by removing the effects of low tSNR voxels in each participant, typically at the disk level, based on the MAP estimate of the tSNR threshold between the spinal cord and CSF/surround voxels (Fig. 3 ). The simultaneous acquisition of brain and spinal cord enabled us to use spatial information from both FOVs for motion correction of the spinal cord EPI volumes. Based on our results reported in Fig. 4 , performing a linear transformation on the entire brain and spinal cord FOVs, followed by slice-wise non-linear correction at each spinal cord slice significantly improves the overall motion correction performance. This is likely due to the 3d linear correction step, which effectively addresses motion correction along the rostro-caudal direction, where the slice-wise correction step is insensitive. Furthermore, FASB registration of the EPI to the template space significantly improves the spinal cord registration compared to current spinal registration approaches (Landelle et al., 2023 ; Vahdat et al., 2015 ). Spinal cord EPI images are often spatially distorted at the disk level, and performing a nonlinear transformation generates non-optimal twisted warping fields. Magnetic susceptibility artifacts also result in signal loss in the distorted voxels, which makes the inclusion of these voxels in further analyses even more problematic. Hence, we proposed performing a simple slice-wise centerline alignment to match the centerlines of the EPI and T2w images. This is possible because we acquire a T2w image at the same angle/slice thickness as the EPI image. Moreover, to standardize the selection of centerline specifically in the distorted slices, we propose to select the EPI centerline based on the estimated tSNR map (Fig. 3 ). Lastly, we account for and remove the effects of distorted slices later by including the signal of participant-specific low tSNR voxels as confound in the group-level regression analysis. Although the T2w medic image is one of the requirements for our registration pipeline, because of its sharp contrast for discerning the gray and white matter, it can also provide useful morphometric information that may be of interest in different clinical populations. For systematic selection of the tSNR threshold between the spinal cord and CSF/surround voxels, we propose to use the optimal boundary (MAP estimate) between two Gaussian distributions, fitted on the tSNR histogram. This makes our threshold selection procedure unbiased and dependent upon the quality of the acquired functional data. In our collected EPI data, the average tSNR threshold was 12.96 (± 0.40 sem). This allowed us to identify the “low tSNR voxels” ( i.e. , those inside the spinal cord mask whose tSNR was lower than the MAP estimate), and remove their effect at the group-level analysis by the inclusion of participant-specific voxel-wise regressors. To the best of our knowledge, our proposed method provides the only available solution to account for and remove the effects of low tSNR voxels from the group-level results, which is particularly critical in spinal cord fMRI due to large magnetic susceptibility artifacts at the spinal cord level. For removal of the scanner noise/physiological artifacts, which are particularly intensified in spinal cord fMRI, we propose to use ICA together with automatic selection of confound components based on their spatial, frequency, and temporal characteristics. Our results show that the performance of ICA and RETROICOR (Glover et al., 2000 ) (based on recorded cardiac and respiratory rates) methods are comparable in terms of reducing the variability attributed to noise (Fig. 6 A-B). Combining the confound regressors extracted using CompCor method (Behzadi et al., 2007 ) (a PCA-based approach) and ICA improves the noise removal performance as evaluated by the explained variance in CSF. However, it should be noted that CSF outcome measures are generally biased toward CompCor , as this method performs PCA inside the CSF mask. Importantly, the analysis of functional connectivity between various spinal cord quadrants (Fig. 6 C) revealed that the various denoising models ( CompCor , ICA, RETROICOR ) did not significantly impact the calculated functional connectivity results. Hence, we included ICA as the method of choice in FASB, due to the ease of acquisition (no need for cardiac and respiratory recordings). However, note that our ICA confound selection procedure is particularly suited for task-based fMRI, in which the frequency and temporal criteria can be defined based on the task timing. In resting-state paradigms selection of confound components can still be achieved (e.g., more active voxels inside the CSF than the spinal cord mask, and larger power outside the neural-related frequency range [0.01 Hz, 0.17 Hz] (Barry et al., 2018a ) than within this range), though it is not as straightforward as in the bask-based design. Hence, for resting-state analysis, a combination of ICA, CompCor or RETROICOR methods may be considered. For validation purposes, we examined the activation maps in the brain and spinal cord using FASB during the performance of a hand force control task. The force control task was selected due to 3 main advantages compared to similar motor tasks performed in the scanner: 1) The task is isometric, so there is no difference in movement kinematics across participants, 2) The produced force is normalized across individuals using their MVC, and 3) The isometric task minimizes motion artifacts, which are critical for spinal cord imaging at the neck level. The activated cortical and subcortical areas during this task are consistent with our previous fMRI studies using hand movements (Braaß et al., 2023 ; Vahdat et al., 2017 , 2015 , 2011 ), including contralateral M1, S1, premotor cortex, SII, PPC, thalamus, striatum, SMA, and ipsilateral cerebellum (Fig. 7 ). Importantly, we found significant task-related activations at the C5-C8 cervical level, mainly at the ipsilateral ventral horn (Fig. 8 ), consistent with the known neuroanatomical location of the finger/wrist motoneurons. Compared to our previous simultaneous spinal cord-brain fMRI studies (Khatibi et al., 2022 ; Kinany et al., 2023 ; Vahdat et al., 2020 , 2015 ), in which inclusion of more than 25 participants were typically recruited to detect significant spinal cord activations, analysis of only 15 participants using the FASB pipeline resulted in significant focal spinal cord activation. Importantly, simultaneous spinal cord-brain fMRI enables investigation of functional connectivity related to the descending and ascending pathways (Landelle et al., 2021 ; Vahdat et al., 2020 ). Our ROI-based analyses demonstrated significant functional connectivity during task performance between cortical, thalamic, cerebellar, and bulbar regions with the cervical ventral horn, and between cortical, thalamic and cerebellar regions with the cervical dorsal horn (Fig. 9 ). Interestingly, we found significant functional connectivity between the ipsilateral medullary reticular formation and the cervical ventral horn, consistent with the involvement of the reticulospinal pathway in hand motor control (Baker, 2011 ; Soteropoulos et al., 2012 ). Since the goal of the current study was to validate FASB pipeline for spinal cord fMRI processing, we performed a simple ROI-based partial correlation analysis to evaluate spinal cord-brain functional connectivity. Partial correlation analysis can account for and remove the effects of task presentation from the BOLD signal time-series, so that their correlation does not simply represent co-activation of areas during task periods as compared to rest (Marrelec et al., 2006 ). Other methods based on ICA (Vahdat et al., 2020 , 2012 ), or BASCO (Rissman et al., 2004 ) (Beta series correlation, particularly suited for an event-related task design) may be used to assess functional connectivity between the brain and spinal cord during task or resting-state conditions. Although our results show that a simple ROI-based correlation analysis may be used to identify significant functional connections in line with the pattern of major descending and ascending projections in humans, future investigations are needed to evaluate the performance of various brain-spinal cord functional connectivity methods. Furthermore, the improved in-plane spatial resolution at the brain level enabled us to detect significant activation clusters in important brainstem nuclei, including the reticular formation (MRF, and PRF), and red nucleus. Importantly, we detected bilateral activation clusters in PRF and MRF, with stronger activation on the ipsilateral side (Fig. 7 ), consistent with the observed pattern of reticulospinal projections targeting forelimb muscles in primates (Sakai et al., 2009 ; Soteropoulos et al., 2012 ). The contribution of various brainstem nuclei in the reticular formation, including the lateral rostral medulla (Arber and Costa, 2022 ; Ruder et al., 2021 ; Ruder and Arber, 2019 ), and caudal medulla (Esposito et al., 2014 ) (corresponding to MRF in humans) in forelimb motor control has been recently the topic of intense investigations in rodents. We hope that FASB acquisition/analysis pipeline provides a standard method for examination of task-related activation and functional interaction across cortical, subcortical, brainstem, and spinal cord levels, and enables testing/translation of preclinical findings to humans and clinical populations. Declarations Funding: Research reported in this study was supported by NIH NINDS under Award Number R01NS131227 (PI: Vahdat), and Florida Department of Health James and Esther King Biomedical Research, Award Number 23K07 (PI: Vahdat). This work was supported in part by an NIH award, S10 OD021726, for High End Instrumentation. JD was supported by the Courtois Foundation, https://www.canadahelps.org/en/charities/fondation-courtois . C.L was supported by Fonds de Recherche du Québec—Santé (FRQ-S) and the Canadian Institutes of Health Research (CIHR). Author Contribution S.V. designed the experiments, optimized the pipeline, and wrote the codes. S.V., F.B., C.L. collected the data, performed data analysis and generated the figures. S.V., O.L., B.D., J.D., F.B., and C.L wrote and revised the manuscript. All authors reviewed the manuscript. References Afyouni, S., Nichols, T.E., 2018. Insight and inference for DVARS. Neuroimage 172, 291–312. https://doi.org/10.1016/J.NEUROIMAGE.2017.12.098 Arber, S., Costa, R.M., 2022. Networking brainstem and basal ganglia circuits for movement. Nat. Rev. Neurosci. 23, 342–360. https://doi.org/10.1038/S41583-022-00581-W Archer, D.B., Kang, N., Misra, G., Marble, S., Patten, C., Coombes, S.A., 2018. Visual feedback alters force control and functional activity in the visuomotor network after stroke. NeuroImage Clin. 17, 505–517. https://doi.org/10.1016/j.nicl.2017.11.012 Archer, D.B., Misra, G., Patten, C., Coombes, S.A., 2016. Microstructural properties of premotor pathways predict visuomotor performance in chronic stroke. Hum. Brain Mapp. 37, 2039–2054. https://doi.org/10.1002/hbm.23155 Baker, S.N., 2011. The primate reticulospinal tract, hand function and functional recovery. J. Physiol. 589, 5603–5612. https://doi.org/10.1113/JPHYSIOL.2011.215160 Barry, R.L., Conrad, B.N., Smith, S.A., Gore, J.C., 2018a. A practical protocol for measurements of spinal cord functional connectivity. Sci. Rep. 8. https://doi.org/10.1038/s41598-018-34841-6 Barry, R.L., Vannesjo, S.J., By, S., Gore, J.C., Smith, S.A., 2018b. Spinal cord MRI at 7T. Neuroimage 168, 437–451. https://doi.org/10.1016/J.NEUROIMAGE.2017.07.003 Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage 37, 90–101. https://doi.org/10.1016/J.NEUROIMAGE.2007.04.042 Braaß, H., Feldheim, J., Chu, Y., Tinnermann, A., Finsterbusch, J., Büchel, C., Schulz, R., Gerloff, C., 2023. Association between activity in the ventral premotor cortex and spinal cord activation during force generation—A combined cortico‐spinal fMRI study. Hum. Brain Mapp. 44, 6471. https://doi.org/10.1002/HBM.26523 Brooks, J.C.W., Beckmann, C.F., Miller, K.L., Wise, R.G., Porro, C.A., Tracey, I., Jenkinson, M., 2008. Physiological noise modelling for spinal functional magnetic resonance imaging studies. Neuroimage 39, 680–692. https://doi.org/10.1016/j.neuroimage.2007.09.018 Burciu, R.G., Chung, J.W., Shukla, P., Ofori, E., Li, H., McFarland, N.R., Okun, M.S., Vaillancourt, D.E., 2016. Functional MRI of disease progression in Parkinson disease and atypical parkinsonian syndromes. Neurology 87, 709–717. https://doi.org/10.1212/WNL.0000000000002985 Cauzzo, S., Singh, K., Stauder, M., García-Gomar, M.G., Vanello, N., Passino, C., Staab, J., Indovina, I., Bianciardi, M., 2022. Functional connectome of brainstem nuclei involved in autonomic, limbic, pain and sensory processing in living humans from 7 Tesla resting state fMRI. Neuroimage 250. https://doi.org/10.1016/J.NEUROIMAGE.2022.118925 Cox, R.W., 1996. AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res. 29, 162–173. https://doi.org/10.1006/cbmr.1996.0014 Davidson, A.G., Schieber, M.H., Buford, J.A., 2007. Bilateral spike-triggered average effects in arm and shoulder muscles from the monkey pontomedullary reticular formation. J. Neurosci. 27, 8053–8058. https://doi.org/10.1523/JNEUROSCI.0040-07.2007 De Leener, B., Fonov, V.S., Collins, D.L., Callot, V., Stikov, N., Cohen-Adad, J., 2018. PAM50: Unbiased multimodal template of the brainstem and spinal cord aligned with the ICBM152 space. Neuroimage 165, 170–179. https://doi.org/10.1016/J.NEUROIMAGE.2017.10.041 De Leener, B., Lévy, S., Dupont, S.M., Fonov, V.S., Stikov, N., Louis Collins, D., Callot, V., Cohen-Adad, J., 2017. SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data. Neuroimage 145, 24–43. https://doi.org/10.1016/j.neuroimage.2016.10.009 Doyon, J., Gabitov, E., Vahdat, S., Lungu, O., Boutin, A., 2018. Current issues related to motor sequence learning in humans. Curr. Opin. Behav. Sci. 20, 89–97. https://doi.org/10.1016/J.COBEHA.2017.11.012 Esposito, M.S., Capelli, P., Arber, S., 2014. Brainstem nucleus MdV mediates skilled forelimb motor tasks. Nature 508, 351–356. https://doi.org/10.1038/nature13023 Finsterbusch, J., Eippert, F., Buchel, C., 2012. Single, slice-specific z-shim gradient pulses improve T2*-weighted imaging of the spinal cord. Neuroimage 59, 2307–2315. https://doi.org/10.1016/j.neuroimage.2011.09.038 Finsterbusch, J., Sprenger, C., Buchel, C., 2013. Combined T2*-weighted measurements of the human brain and cervical spinal cord with a dynamic shim update. Neuroimage 79, 153–161. https://doi.org/10.1016/j.neuroimage.2013.04.021 Glover, G.H., Li, T.Q., Ress, D., 2000. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn. Reson. Med. 44, 162–167. https://doi.org/10.1002/1522-2594(200007)44:13.0.CO;2-E Gros, C., De Leener, B., Badji, A., Maranzano, J., Eden, D., Dupont, S.M., Talbott, J., Zhuoquiong, R., Liu, Y., Granberg, T., Ouellette, R., Tachibana, Y., Hori, M., Kamiya, K., Chougar, L., Stawiarz, L., Hillert, J., Bannier, E., Kerbrat, A., Edan, G., Labauge, P., Callot, V., Pelletier, J., Audoin, B., Rasoanandrianina, H., Brisset, J.C., Valsasina, P., Rocca, M.A., Filippi, M., Bakshi, R., Tauhid, S., Prados, F., Yiannakas, M., Kearney, H., Ciccarelli, O., Smith, S., Treaba, C.A., Mainero, C., Lefeuvre, J., Reich, D.S., Nair, G., Auclair, V., McLaren, D.G., Martin, A.R., Fehlings, M.G., Vahdat, S., Khatibi, A., Doyon, J., Shepherd, T., Charlson, E., Narayanan, S., Cohen-Adad, J., 2019. Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks. Neuroimage 184, 901–915. https://doi.org/10.1016/j.neuroimage.2018.09.081 Hyvarinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Netw 10, 626–634. https://doi.org/10.1109/72.761722 Islam, H., Law, C.S.W., Weber, K.A., Mackey, S.C., Glover, G.H., 2019. Dynamic per slice shimming for simultaneous brain and spinal cord fMRI. Magn. Reson. Med. 81, 825–838. https://doi.org/10.1002/MRM.27388 Kaptan, M., Horn, U., Vannesjo, S.J., Mildner, T., Weiskopf, N., Finsterbusch, J., Brooks, J.C.W., Eippert, F., 2023. Reliability of resting-state functional connectivity in the human spinal cord: Assessing the impact of distinct noise sources. Neuroimage 275. https://doi.org/10.1016/J.NEUROIMAGE.2023.120152 Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. J. Neurosci. Methods 276, 56–72. https://doi.org/10.1016/J.JNEUMETH.2016.10.019 Khatibi, A., Vahdat, S., Lungu, O., Finsterbusch, J., Büchel, C., Cohen-Adad, J., Marchand-Pauvert, V., Doyon, J., 2022. Brain-spinal cord interaction in long-term motor sequence learning in human: An fMRI study. Neuroimage 253, 119111. https://doi.org/10.1016/J.NEUROIMAGE.2022.119111 Kinany, N., Khatibi, A., Lungu, O., Finsterbusch, J., Büchel, C., Marchand-Pauvert, V., Van De Ville, D., Vahdat, S., Doyon, J., 2023. Decoding cerebro-spinal signatures of human behavior: Application to motor sequence learning. Neuroimage 275. https://doi.org/10.1016/J.NEUROIMAGE.2023.120174 Kinany, N., Pirondini, E., Micera, S., Van De Ville, D., 2022. Spinal Cord fMRI: A New Window into the Central Nervous System. Neuroscientist. https://doi.org/10.1177/10738584221101827 Kong, Y., Jenkinson, M., Andersson, J., Tracey, I., Brooks, J.C.W., 2012. Assessment of physiological noise modelling methods for functional imaging of the spinal cord. Neuroimage 60, 1538–1549. https://doi.org/10.1016/j.neuroimage.2011.11.077 Landelle, C., Dahlberg, L.S., Lungu, O., Misic, B., De Leener, B., Doyon, J., 2023. Altered Spinal Cord Functional Connectivity Associated with Parkinson’s Disease Progression. Mov. Disord. 38, 636–645. https://doi.org/10.1002/MDS.29354 Landelle, C., Lungu, O., Vahdat, S., Kavounoudias, A., Marchand-Pauvert, V., De Leener, B., Doyon, J., 2021. Investigating the human spinal sensorimotor pathways through functional magnetic resonance imaging. Neuroimage 245. https://doi.org/10.1016/J.NEUROIMAGE.2021.118684 Lemon, R.N., 2008. Descending pathways in motor control. Annu Rev Neurosci 31, 195–218. https://doi.org/10.1146/annurev.neuro.31.060407.125547 Marrelec, G., Krainik, A., Duffau, H., Pelegrini-Issac, M., Lehericy, S., Doyon, J., Benali, H., 2006. Partial correlation for functional brain interactivity investigation in functional MRI. Neuroimage 32, 228–237. https://doi.org/10.1016/j.neuroimage.2005.12.057 Pierrot-Deseilligny, E., Burke, D., 2012. The circuitry of the human spinal cord: Spinal and corticospinal mechanisms of movement, The Circuitry of the Human Spinal Cord: Spinal and Corticospinal Mechanisms of Movement. https://doi.org/10.1017/CBO9781139026727 Raudino, F., Leva, S., 2012. Involvement of the spinal cord in Parkinson’s disease. Int. J. Neurosci. 122, 1–8. https://doi.org/10.3109/00207454.2011.613551 Rissman, J., Gazzaley, A., D’Esposito, M., 2004. Measuring functional connectivity during distinct stages of a cognitive task. Neuroimage 23, 752–763. https://doi.org/10.1016/J.NEUROIMAGE.2004.06.035 Ruder, L., Arber, S., 2019. Brainstem Circuits Controlling Action Diversification. Annu. Rev. Neurosci. 42, 485–504. https://doi.org/10.1146/ANNUREV-NEURO-070918-050201 Ruder, L., Schina, R., Kanodia, H., Valencia-Garcia, S., Pivetta, C., Arber, S., 2021. A functional map for diverse forelimb actions within brainstem circuitry. Nature 590, 445–450. https://doi.org/10.1038/S41586-020-03080-Z Sakai, S.T., Davidson, A.G., Buford, J.A., 2009. Reticulospinal neurons in the pontomedullary reticular formation of the monkey (Macaca fascicularis). Neuroscience 163, 1158–1170. https://doi.org/10.1016/J.NEUROSCIENCE.2009.07.036 Singh, K., García-Gomar, M.G., Bianciardi, M., 2021. Probabilistic Atlas of the Mesencephalic Reticular Formation, Isthmic Reticular Formation, Microcellular Tegmental Nucleus, Ventral Tegmental Area Nucleus Complex, and Caudal-Rostral Linear Raphe Nucleus Complex in Living Humans from 7 Tesla Magnetic Resonance Imaging. Brain Connect. 11, 613–623. https://doi.org/10.1089/BRAIN.2020.0975 Smith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E., Johansen-Berg, H., Bannister, P.R., De Luca, M., Drobnjak, I., Flitney, D.E., Niazy, R.K., Saunders, J., Vickers, J., Zhang, Y., De Stefano, N., Brady, J.M., Matthews, P.M., 2004. Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23 Suppl 1, S208-19. https://doi.org/S1053-8119(04)00393-3 [pii]10.1016/j.neuroimage.2004.07.051 Soteropoulos, D.S., Williams, E.R., Baker, S.N., 2012. Cells in the monkey ponto-medullary reticular formation modulate their activity with slow finger movements. J. Physiol. 590, 4011–4027. https://doi.org/10.1113/jphysiol.2011.225169 Tinnermann, A., Büchel, C., Cohen-Adad, J., 2021. Cortico-spinal imaging to study pain. Neuroimage. https://doi.org/10.1016/j.neuroimage.2020.117439 Vahdat, S., Darainy, M., Milner, T.E., Ostry, D.J., 2011. Functionally specific changes in resting-state sensorimotor networks after motor learning. J Neurosci 31, 16907–16915. https://doi.org/10.1523/JNEUROSCI.2737-11.2011 Vahdat, S., Darainy, M., Ostry, D.J., 2014. Structure of Plasticity in Human Sensory and Motor Networks Due to Perceptual Learning. J. Neurosci. 34, 2451–2463. https://doi.org/10.1523/JNEUROSCI.4291-13.2014 Vahdat, S., Fogel, S., Benali, H., Doyon, J., 2017. Network-wide reorganization of procedural memory during NREM sleep revealed by fMRI. Elife 6, e24987. https://doi.org/10.7554/eLife.24987 Vahdat, S., Khatibi, A., Lungu, O., Finsterbusch, J., Büchel, C., Cohen-Adad, J., Marchand-Pauvert, V., Doyon, J., 2020. Resting-state brain and spinal cord networks in humans are functionally integrated. PLoS Biol. 18(7): e30. https://doi.org/10.1371/journal.pbio.3000789 Vahdat, S., Lungu, O., Cohen-Adad, J., Marchand-Pauvert, V., Benali, H., Doyon, J., 2015. Simultaneous brain–cervical cord fMRI reveals intrinsic spinal cord plasticity during motor sequence learning. PLoS Biol. 13, 1–25. https://doi.org/10.1371/journal.pbio.1002186 Vahdat, S., Maneshi, M., Grova, C., Gotman, J., Milner, T.E., 2012. Shared and specific independent components analysis for between-group comparison. Neural Comput. 24, 3052–3090. https://doi.org/10.1162/NECO_a_00355 Vaillancourt, D.E., Thulborn, K.R., Corcos, D.M., 2003. Neural Basis for the Processes That Underlie Visually Guided and Internally Guided Force Control in Humans. J. Neurophysiol. 90, 3330–3340. https://doi.org/10.1152/jn.00394.2003 Wahl, A.S., Schwab, M.E., 2014. Finding an optimal rehabilitation paradigm after stroke: Enhancing fiber growth and training of the brain at the right moment. Front. Hum. Neurosci. 8, 381. https://doi.org/10.3389/FNHUM.2014.00381/BIBTEX Additional Declarations No competing interests reported. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3889284","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":269902678,"identity":"17d60d8e-4d7a-40bd-b552-d7a420df5982","order_by":0,"name":"Shahabeddin Vahdat","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAtElEQVRIiWNgGAWjYHACNgYeBgY5BubDDSA20VoMjBnYEknUkthAtBb5GcnPHryp+ZO+4RhjA8OHssOEtRjcSDM3nHPMIBekhXHGOWK08Bwwk+ZhA2q539jAzNtGhBb5nuPfpHn+GaQbAG1h/kuMFobjPWbSvG0GCWAtjMRoMTjeUyY5t8/YcCZQy8Gec+lEOKyZfZvEm29y8nzHmA8++FFmTYTDkMEBEtWPglEwCkbBKMAFAGEuOn4Sa8NyAAAAAElFTkSuQmCC","orcid":"","institution":"University of Florida","correspondingAuthor":true,"prefix":"","firstName":"Shahabeddin","middleName":"","lastName":"Vahdat","suffix":""},{"id":269902683,"identity":"f288a0e8-ef48-4c50-93cd-1c48f05c4aca","order_by":1,"name":"Caroline Landelle","email":"","orcid":"","institution":"McGill University","correspondingAuthor":false,"prefix":"","firstName":"Caroline","middleName":"","lastName":"Landelle","suffix":""},{"id":269902686,"identity":"f8d4afeb-a04d-4c8d-afbe-a9b27c200aac","order_by":2,"name":"Ovidiu Lungu","email":"","orcid":"","institution":"McGill University","correspondingAuthor":false,"prefix":"","firstName":"Ovidiu","middleName":"","lastName":"Lungu","suffix":""},{"id":269902687,"identity":"a47a46ba-76d0-4fe7-a13b-0b330091134b","order_by":3,"name":"Benjamin De Leener","email":"","orcid":"","institution":"Polytechnique Montreal","correspondingAuthor":false,"prefix":"","firstName":"Benjamin","middleName":"","lastName":"De Leener","suffix":""},{"id":269902688,"identity":"2f732d67-e7d5-45f7-9b94-50eb754b259b","order_by":4,"name":"Julien Doyon","email":"","orcid":"","institution":"McGill University","correspondingAuthor":false,"prefix":"","firstName":"Julien","middleName":"","lastName":"Doyon","suffix":""},{"id":269902690,"identity":"56c7bbe4-ab9a-49bd-9b91-4f3db8b36b35","order_by":5,"name":"Fatemeh Baniasad","email":"","orcid":"","institution":"University of Florida","correspondingAuthor":false,"prefix":"","firstName":"Fatemeh","middleName":"","lastName":"Baniasad","suffix":""}],"badges":[],"createdAt":"2024-01-22 23:44:05","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3889284/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3889284/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":50390194,"identity":"fd40671d-4267-4d1e-9b08-9c3cff35645f","added_by":"auto","created_at":"2024-01-30 18:40:13","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":368829,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFASB EPI analysis pipeline.\u003c/strong\u003e The main steps for analysis of EPI spinal cord data include: 2-step motion correction, GRE field map correction, slice-wise centerline alignment to T2w image using tSNR map as a guide, MAP estimation of tSNR threshold, spatial smoothing, noise removal by regressing out motion correction parameters, average CSF signal, and ICA-based confounds, high-pass temporal filtering, individual-level regression using task timing as the regressor of interest and motion outliers, registration of the summary statistics images to the spinal cord template, group-level regression by inclusion of individual-specific voxel-wise low-tSNR confounds, and correction for multiple comparisons using GRF.\u003c/p\u003e","description":"","filename":"floatimage1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/24bd43f39d88258a750928d4.jpg"},{"id":50389820,"identity":"2368c6c5-1c80-40cf-815a-c5d573d44a7b","added_by":"auto","created_at":"2024-01-30 18:32:13","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":314505,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFASB registration pipeline.\u003c/strong\u003e EPI spinal cord images are registered to the spinal cord template (PAM50) in 3 steps: 1) the reference EPI volume is registered to the T2w image using slice-wise centerline alignment, 2) T2w image is registered to T1w image via a slice-wise rigid body transformation using \u003cem\u003esct_register_multimodal \u003c/em\u003e(SCT), and 3) T1w image is registered to PAM50 via a nonlinear slice-wise transformation using \u003cem\u003esct_register_to_template\u003c/em\u003e (SCT).\u003c/p\u003e","description":"","filename":"floatimage2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/e4d906042326283c437c9241.jpg"},{"id":50389209,"identity":"77f3edeb-c4d9-4bc3-ba2a-b80e98265faa","added_by":"auto","created_at":"2024-01-30 18:24:13","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":307315,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe tSNR threshold estimation procedure. \u003c/strong\u003eTop left: the tSNR map in the spinal cord and the surround mask for a representative participant. Some slices around the disk level include low tSNR voxels in the spinal cord due to exacerbated susceptibility artifacts. Top right: the average tSNR values across all participants inside the spinal cord mask, the surround mask (containing mainly CSF), and the spinal cord tSNR-corrected mask. Bottom: the distribution of the tSNR values inside the spinal cord and surround mask for a representative participant, fitted using a Gaussian mixture distribution model using two components. The dashed line represents the \u003cem\u003etSNR threshold,\u003c/em\u003e defined as the Maximum a Posteriori (MAP) estimate of the optimal boundary between the two Gaussian distributions. tSNR: temporal Signal-Noise Ratio\u003c/p\u003e","description":"","filename":"floatimage3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/66ab1d767f2326aed91413d0.jpg"},{"id":50389821,"identity":"dbd82b5c-b8ca-4e76-b32a-a908065b21d7","added_by":"auto","created_at":"2024-01-30 18:32:13","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":102883,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePerformance of various motion correction procedures for spinal cord EPI data.\u003c/strong\u003eLeft and right panels, respectively, show the average performance across participants of different motion correction procedures as evaluated by DVARS and tSNR in the spinal cord mask. The combined linear and non-linear procedure (FASB) outperforms the other methods [linear (SP+BR Flirt), and nonlinear (Slice-wise SCT)] in terms of both DVARS and tSNR. “Raw” reports results without motion correction. * and ** indicates \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.05 and \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.01, repeated-measures ANOVA. Error bar: sem. tSNR: temporal Signal-Noise Ratio, SP: Spinal cord; BR: Brain\u003c/p\u003e","description":"","filename":"floatimage4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/1359cebb4002caa42a026068.jpg"},{"id":50389208,"identity":"d6e1ff42-9c53-44c1-927c-d969bba90677","added_by":"auto","created_at":"2024-01-30 18:24:13","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":137035,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePerformance of various registration methods to align the EPI spinal cord image to the template. \u003c/strong\u003eLeft and right panels, respectively, show the average performance across participants of different registration methods as evaluated by Error and Overlap ratios associated with the registered spinal cord mask in the PAM50 space. FASB registration method outperforms the other methods [T1w-based registration, and GRE correction+T1w-based registration] in terms of both Error and Overlap measures. * and ** indicates \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.05 and \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.001, repeated-measures ANOVA. Error bar: sem.\u003c/p\u003e","description":"","filename":"floatimage5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/d99492f569f9fc45cf6daf5d.jpg"},{"id":50389211,"identity":"b5dc5221-5f1e-4737-a758-06c37d6cb35c","added_by":"auto","created_at":"2024-01-30 18:24:13","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":222051,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePerformance of various noise removal models on the collected resting-state fMRI data.\u003c/strong\u003e \u003cstrong\u003e(A)\u003c/strong\u003e shows the ratio between the R\u003csup\u003e2\u003c/sup\u003e of CSF and spinal cord using different noise removal models, averaged across participants. \u003cem\u003eCompCor\u003c/em\u003e outperforms the other models in terms of R\u003csup\u003e2\u003c/sup\u003e ratio. \u003cstrong\u003e(B)\u003c/strong\u003e shows the adjusted R\u003csup\u003e2\u003c/sup\u003e (AdjR\u003csup\u003e2\u003c/sup\u003e) inside the CSF mask, averaged across participants. \u003cem\u003eICA\u003c/em\u003e, \u003cem\u003eCompCor\u003c/em\u003e and \u003cem\u003eICA + CompCor\u003c/em\u003e models outperform the \u003cem\u003emeanCSF\u003c/em\u003e and \u003cem\u003eRetroicor\u003c/em\u003e models in terms of AdjR\u003csup\u003e2\u003c/sup\u003e. \u003cstrong\u003e(C)\u003c/strong\u003e shows functional connectivity between various spinal cord quadrants using various denoising models. As shown, the choice of model does not significantly impact the calculated functional connectivity results. ** indicates \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.001. Error bar: sem.\u003c/p\u003e","description":"","filename":"floatimage6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/a8177615cf7300cd56280547.jpg"},{"id":50389216,"identity":"7c17c526-54e1-4374-8bc8-17b2af8b21bf","added_by":"auto","created_at":"2024-01-30 18:24:13","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":356379,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eBrain activation map (in radiological view) associated with unilateral handgrip force control.\u003c/strong\u003e \u003cstrong\u003e(A)\u003c/strong\u003e Unilateral handgrip force control (right hand) results in significant activation clusters mainly in the contralateral sensorimotor cortical, thalamic, and striatal areas, as well as ipsilateral cerebellar and brainstem regions. \u003cstrong\u003e(B)\u003c/strong\u003e The zoomed-in views of sagittal slices covering the brainstem highlight significant activation clusters in the ipsilateral MRF and PRF regions. The anatomical locations of MRF and PRF are respectively highlighted in blue and green colors. PMd: dorsal premotor cortex, SMA: supplementary motor area, M1: primary motor cortex, S1: primary somatosensory cortex, SII: secondary somatosensory cortex, PUT: putamen, THL: thalamus, RN: red nucleus, CB IV-VI: cerebellum lobule IV-VI, PRF: pontine reticular formation, MRF: medullary reticular formation. R: right (ipsilateral), L: left (contralateral). All activation maps are corrected for multiple comparisons using GRF, \u003cem\u003eZ\u003c/em\u003e \u0026gt; 2.7, \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.05.\u003c/p\u003e","description":"","filename":"floatimage7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/d2403ab688b799be283d58f3.jpg"},{"id":50389823,"identity":"9fe8f2d6-b423-4972-a016-ec259eda2231","added_by":"auto","created_at":"2024-01-30 18:32:13","extension":"jpg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":235569,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSpinal cord activation map associated with unilateral handgrip force control using various analysis methods.\u003c/strong\u003e From left to right, each panel shows the significant activation clusters associated with unilateral handgrip force control using: FASB (ICA inclusion and tSNR correction), only ICA inclusion (ICA-NotSNR), only tSNR correction (NoICA-tSNR), and no ICA inclusion nor tSNR correction (NoICA-NotSNR). As expected, the activation clusters are located in the ipsilateral ventral horn at the lower cervical levels. The FASB method results in the greatest number of activation clusters covering the C5-C8 cervical segmental levels. All activation maps are corrected for multiple comparisons using GRF, \u003cem\u003eZ\u003c/em\u003e \u0026gt; 2.7, \u003cem\u003ep\u003c/em\u003e \u0026lt; 0.05. tSNR: temporal Signal-Noise Ratio, ICA: Independent Component Analysis.\u003c/p\u003e","description":"","filename":"floatimage8.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/23a0012a2170c31d4d779ef1.jpg"},{"id":50389214,"identity":"df36d0c5-42d5-4398-a0cd-68e68bc576e2","added_by":"auto","created_at":"2024-01-30 18:24:13","extension":"jpg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":154179,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional connectivity between various brain and spinal cord regions calculated using FASB pipeline.\u003c/strong\u003e Figure shows significant functional connectivity of various brain regions with the ipsilateral ventral and dorsal cervical cord (C5-C8 levels) as evaluated by partial correlation analysis during unilateral handgrip force control task. Postcent: postcentral gyrus, Precent: precentral gyrus, Thal: thalamus, CB V-VI: cerebellum lobule V-VI, SMA: supplementary motor area, MRF: medullary reticular formation, VenH: ventral horn, DorH: dorsal horn, R: right (ipsilateral), L: left (contralateral). \u003cem\u003ep-\u003c/em\u003evalues reports repeated measures \u003cem\u003et\u003c/em\u003e-statistics.\u003c/p\u003e","description":"","filename":"floatimage9.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/15d34e1dc2f62db661c1c5ce.jpg"},{"id":51364830,"identity":"744e9426-e8f4-466f-ae56-1532cd35d53d","added_by":"auto","created_at":"2024-02-20 09:45:30","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1384274,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3889284/v1/ee737762-bce3-4e4c-bf07-8818076b575a.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"FASB: an integrated processing pipeline for Functional Analysis of simultaneous Spinal cord-Brain fMRI","fulltext":[{"header":"Introduction","content":"\u003cp\u003eA plethora of evidence suggests that many of the movement disorders, such as amyotrophic lateral sclerosis, multiple sclerosis, Parkinson\u0026rsquo;s disease, stroke, spinal cord injury, dystonia, and cerebral palsy impair both the brain and spinal cord circuits (Pierrot-Deseilligny and Burke, \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Raudino and Leva, \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Wahl and Schwab, \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). However, human neuroimaging studies with patients affected by these disorders have mainly focused on characterizing the brain circuit dysfunction and deficits, while disregarding the critical role of the spinal cord circuitry and brain-spinal cord pathways. Simultaneous functional magnetic resonance imaging (fMRI) of the spinal cord and brain provides an unprecedented opportunity to address this limitation by allowing the assessment of the ascending sensory and descending motor pathways functioning in humans \u003cem\u003ein vivo\u003c/em\u003e (Kinany et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Tinnermann et al., \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). By employing this combined approach, researchers can investigate not only the neural activations at multiple levels of the central nervous system (CNS), including cortical, subcortical, brainstem, and spinal cord, but also functional interactions among these structures, thereby gaining valuable insights into the neural mechanisms underlying motor and sensory processing. Nevertheless, the imaging parameters and processing pipeline associated with this simultaneous neuroimaging approach are still in the early stages of development. In this study, we introduce an optimized acquisition protocol including the required anatomical and functional images, as well as an integrated analysis pipeline, named \u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eF\u003c/span\u003eunctional \u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eA\u003c/span\u003enalysis of simultaneous \u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eS\u003c/span\u003epinal cord-\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eB\u003c/span\u003erain fMRI (FASB). The proposed pipeline aims to improve the current challenges associated with the processing of spinal cord fMRI data under simultaneous spinal cord-brain acquisition in three domains. These improvements include motion correction, registration to the template, and removing the effects of large susceptibility artifacts, often exacerbated at the spinal disk levels.\u003c/p\u003e \u003cp\u003ePrior work using simultaneous spinal cord-brain fMRI has revealed that task- and context-dependent dynamical interactions between brain regions and the spinal cord subserve various sensory and motor functions. In the last decade, the optimization of acquisition protocols at high-field strength (Barry et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2018b\u003c/span\u003e), as well as the development of multi-slab pulse sequences designed to improve the acquisition of selective fields-of-view (FOVs), combined with advanced shimming procedures (Finsterbusch et al., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2013\u003c/span\u003e, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Islam et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2019\u003c/span\u003e) have significantly enhanced the quality of spinal cord-brain fMRI. Using these advanced image acquisition and shimming techniques, we have recently demonstrated the existence of integrated spinal cord and brain resting-state networks that follow neuroanatomical principles governing the ascending and descending pathways in humans (Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Despite clear advancement in spinal cord fMRI thanks to these new techniques, there is still work needed to improve image quality due to all the challenges associated with functional imaging of the spinal cord (Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eThe application of spinal cord fMRI is starting to augment our brain-centered view of neurological disorders, by investigating the role of spinal cord circuit function in movement disorders, \u003cem\u003ee.g.\u003c/em\u003e in Parkinson\u0026rsquo;s disease (Landelle et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Yet studying the function of ascending and descending pathways requires the application of simultaneous fMRI methods, which has been limited in clinical populations due to a number of technical challenges (Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). These include the complex acquisition procedures that can achieve both a reasonable signal to noise ratio (SNR), and a high spatial resolution for distinguishing spinal cord quadrants, further constrained by the scanning time considerations. The goal of this study is to propose a simple acquisition method that can be implemented in clinical populations, together with an integrated processing pipeline to standardize the procedures for simultaneous spinal cord-brain fMRI.\u003c/p\u003e \u003cp\u003eThe spinal cord toolbox (SCT) (De Leener et al., \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2017\u003c/span\u003e) has offered optimized set of tools for preprocessing of the spinal cord MRI data, including slice-wise motion correction, automatic segmentation of the spinal cord using convolutional neural networks (Gros et al., \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), and nonlinear slice-wise registration to the template (De Leener et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). However, the combined acquisition of the spinal cord and brain provides new opportunities, as well as challenges, which can benefit from a customized processing pipeline. For example, the motion correction procedure implemented in the SCT is solely based on the spinal cord FOV, and can thus be improved by incorporating information from the acquired brain FOV. Also, large spatial distortions due to exacerbated susceptibility artifacts at the disk levels in BOLD imaging limit the efficacy of non-linear transformations for registration. Importantly, this may cause signal loss in the EPI data and may degrade SNR in some of the spinal cord slices. These low tSNR voxels practically degrade the group-level results by adding noisy voxels to the group-level averaging. Yet, there are two solutions to the latter problem. One approach is to remove the low SNR voxels from the participant\u0026rsquo;s spinal cord mask. However, the voxels included at the group-level analysis should be present in the mask of every single individuals. Because low-SNR voxels appear at different locations across participants (depending on neck curvature, surrounding tissues, etc.), a large number of voxels will be removed from the group-level analysis, which makes this approach impractical. Another approach is to remove the effect of these low-SNR voxels statistically from the group-level analysis, by incorporating participant-specific voxel-wise confounds in the group-level regression model. In this study, we employ the latter idea and propose a new approach that features automated modeling and identification of spinal cord voxels with low temporal SNR in each participant, followed by removing their effects statistically at the group-level analysis without the need to exclude those voxels from the entire group.\u003c/p\u003e \u003cp\u003eTo evaluate the efficacy of the proposed pipeline, FASB, we collected simultaneous spinal cord-brain fMRI data during the performance of a motor control task (unilateral handgrip force control), as well as during resting-state conditions. First, we compared the performance of the proposed algorithm for motion correction, registration to the spinal cord template, and physiological noise removal, to the current methods used in the spinal cord fMRI analysis. We then reported the activation maps in the brain, brainstem and spinal cord using FASB and other analysis methods during the motor control task. We hypothesized that the subject\u0026rsquo;s performance of this unilateral handgrip task would significantly increase the level of activity in the ipsilateral ventral horn at the C6-C8 cervical levels (Pierrot-Deseilligny and Burke, \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2012\u003c/span\u003e), as well as in the contralateral primary and secondary sensorimotor cortical areas, striatum, thalamus, ipsilateral cerebellar cortex (Doyon et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2011\u003c/span\u003e), and ipsilateral motor-related brainstem nuclei, including pontine and medullary reticular formation (Arber and Costa, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Ruder and Arber, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). We further assessed the functional connectivity of the cervical dorsal and ventral horns with various sensory and motor areas of the brain and brainstem. Our results demonstrate that the proposed acquisition protocol and processing pipeline (FASB) outperforms other methods commonly used in spinal cord fMRI analysis, and provides a robust method for accurately characterizing the activation and functional connectivity of distinct cortical, subcortical, brainstem and spinal cord regions in humans in vivo.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eParticipants\u003c/h2\u003e \u003cp\u003e This study was reviewed and approved by the ethics committee at the University of Florida, which adhered to the Declaration of Helsinki. All participants signed an informed consent form prior to participating in the study and were debriefed and compensated at the end of the experiment. Fifteen young, right-handed, healthy adults (10 females, mean age\u0026thinsp;=\u0026thinsp;22.7 years, SD\u0026thinsp;=\u0026thinsp;4.0) were selected to take part in this study based on the following exclusion criteria: presence of a history of neurological and psychiatric diseases, of any motor-system complication, the current use of any medication, and presence of any MRI-incompatible object in the body. An additional twenty-one healthy participants were recruited for resting-state acquisition (12 females, mean age\u0026thinsp;=\u0026thinsp;31.4 years, SD\u0026thinsp;=\u0026thinsp;7.1). These participants fulfilled the inclusion criteria previously described and the data acquisition was approved by the ethics committee at the University of Montr\u0026eacute;al, which also adheres to the Declaration of Helsinki. Due to poor fMRI image quality, resting-state data from two participants from this sample were not included in analyses, and data from another participant was also excluded due to missing physiological noise measurements. Thus, the resting-state sample included in the final analyses consisted of 18 healthy participants (mean age 32.39 (SD\u0026thinsp;=\u0026thinsp;13.19), F/M\u0026thinsp;=\u0026thinsp;11/7).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eMotor control task\u003c/h2\u003e \u003cp\u003eParticipants who performed a motor control task in the scanner were asked to produce isometric hand grip force with their right hand using an MR-compatible fiber-optic force transducer. Prior to the MRI scan, participants were first trained on the task, and the grip maximum voluntary contraction (MVC) was measured. Specifically, they were asked to produce unilateral grip force using their right hand by opposing their thumb, index, and middle fingers against the force transducer (Archer et al., \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2018\u003c/span\u003e, \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Burciu et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). Two runs were completed in a block-design paradigm. Each run consisted of 9 blocks of 15 sec force control, followed by 25 sec rest periods (6 min and 25 sec in total). Throughout the run, 2 bars were displayed on a video monitor that participants saw through a mirror mounted on the head coil: a target bar and a force bar. A change in color of the force bar was used to cue participants to either push or release the force sensor. Green was a go signal for producing and sustaining force, while red indicated a rest period. The target bar was set at 10, 20, or 30% of MVC for each participant (3 blocks at each level, pseudo-randomly ordered), while the force bar moved in the vertical plane according to the force output. During the go signal, participants were required to adjust their grip force levels to match the force and target bars. Force was measured using custom-designed load cells (Vaillancourt et al., \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2003\u003c/span\u003e) and a Micron Optics fiber optic interrogator.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eMRI data acquisition\u003c/h2\u003e \u003cp\u003e \u003cb\u003eMotor control task dataset\u003c/b\u003e. A 3T Siemens Prisma scanner equipped with a 64-channel head/neck coil was used for imaging subjects. We acquired functional magnetic resonance imaging (fMRI) data, covering the whole brain and cervical spinal cord (first cervical, C1, to first thoracic, T1, vertebrae), using blood oxygenation level dependent (BOLD) contrast via a simultaneous multi-slice (SMS) echo-planar imaging (EPI) sequence. Participants laid on the scanner table in a supine position with their head and neck fully supported using foam pads to minimize their motion. We used the neck and brachial plexus SatPad\u0026trade; sets (SatPad, clinical imaging solutions) around the neck to minimize inhomogeneity artifacts and enhance spinal cord functional image quality. Based on previous work in the lab, we recommend their use in fMRI studies that target the spinal cord. Participants were placed in the scanner such that the mid-chin level was in the scanner\u0026rsquo;s isocenter, corresponding to C2-C3 vertebral level on the axial plane. The following parameters were used for EPI acquisitions: 72 axial slices angled orthogonal to the cord at C5 vertebral level, in-plane resolution: 1.6x1.6 mm\u003csup\u003e2\u003c/sup\u003e, matrix size: 120x120, 50% phase oversampling, slice thickness: 4 mm (no gap), TR: 1.85 s, 210 volumes, TE: 27 ms, multiband acceleration factor: 3, GRAPPA factor: 2, 0.875 partial Fourier encoding, FA: 70˚, PE direction: A-P. Additionally, 2 saturation pulses were applied in a V-shaped configuration to minimize ghosting and inflow artifacts related to blood flow in the major cervical vessels.\u003c/p\u003e \u003cp\u003eA 3D-MPRAGE sequence was used to acquire T1-weighted anatomical images covering the entire head as well as neck down to the T3 vertebral level using the following parameters: field of view (FOV): 192\u0026times;240\u0026times;320 mm\u003csup\u003e3\u003c/sup\u003e, sagittal slices; TR: 2300 ms, TE: 3.41 ms, and resolution: 1\u0026times;1\u0026times;1 mm\u003csup\u003e3\u003c/sup\u003e. In addition, a multi-echo data image combination (MEDIC) sequence was used to acquire a T2*-weighted (T2w) anatomical image using the following parameters: axial slices\u0026thinsp;=\u0026thinsp;36, slice thickness\u0026thinsp;=\u0026thinsp;4 mm (no gap), in-plane resolution: 0.75\u0026times;0.75 mm\u003csup\u003e2\u003c/sup\u003e, matrix size: 256\u0026times;256, TR: 1280 ms, TE: 18, 20, 21, and 22 ms; FA: 20˚, and GRAPPA acceleration factor: 2. Although the orientation and the size of slices were matched between the T2w and EPI images, the T2w image had fewer slices, covering only the cervical cord to keep the acquisition time down to 12 min. Importantly, the position of the T2w image FOV in the dorsocaudal orientation was set such that the position of the most caudal slice matches that of the EPI image.\u003c/p\u003e \u003cp\u003e\u003cb\u003eResting-state dataset.\u003c/b\u003e Data were also acquired using a 3-Tesla MRI Scanner (Magnetom-Prisma, Siemens, Erlangen, Germany) equipped with a 64-channel head (1\u0026ndash;7 elements active) and neck coil (1\u0026ndash;2 elements active). T2-weighted images were collected in the transverse orientation with the following parameters: TR\u0026thinsp;=\u0026thinsp;33ms; TE\u0026thinsp;=\u0026thinsp;14ms; FOV\u0026thinsp;=\u0026thinsp;224 mm x 258 mm; flip angle\u0026thinsp;=\u0026thinsp;5\u003csup\u003eo\u003c/sup\u003e; in-plane voxel resolution\u0026thinsp;=\u0026thinsp;0.35 x 0.35 mm\u003csup\u003e2\u003c/sup\u003e, slice thickness\u0026thinsp;=\u0026thinsp;2 mm. A total of 112 slices were acquired, spanning from the top of the cerebellum to the upper thoracic region (approximately at T1), thereby encompassing the entirety of the cervical spinal cord. Resting-state functional images were acquired using an echo-planar imaging (EPI) gradient echo sequence covering the brain and cervical spinal (TR/TE\u0026thinsp;=\u0026thinsp;1550/23 ms, axial FOV\u0026thinsp;=\u0026thinsp;120 x 120 mm\u003csup\u003e2\u003c/sup\u003e, flip angle\u0026thinsp;=\u0026thinsp;70\u0026deg;, in-plane voxel resolution\u0026thinsp;=\u0026thinsp;1.6 x 1.6 mm\u003csup\u003e2\u003c/sup\u003e; slice thickness\u0026thinsp;=\u0026thinsp;4 mm (no gap); iPAT acceleration factor PE\u0026thinsp;=\u0026thinsp;2, iPAT acceleration factor slice\u0026thinsp;=\u0026thinsp;3). A total of 69 axial slices were acquired covering the top of the head to approximately T1 vertebra. Only the slices containing the cervical spinal cord were included in further analyses. One resting state-run (\u003cem\u003ei.e.\u003c/em\u003e, no explicit task) was acquired in total with an acquisition time of 10 minutes and 35 seconds each (230 volumes). Participants were instructed to relax watch the video and minimize swallowing and motion. Physiological recordings were acquired using a pulse sensor and a respiration belt (Siemens Physiology Monitoring Unit).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eSpinal cord analysis pipeline\u003c/h2\u003e \u003cp\u003eWe utilized bash programming, the FMRIB Software Library (FSL) (Smith et al., \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2004\u003c/span\u003e), Spinal Cord Toolbox (SCT version 6.0) (De Leener et al., \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2017\u003c/span\u003e), and in-house MATLAB programs to process the collected data. The full spinal cord processing pipeline will be publicly available at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e\u003ca href=\"http://www.vahdatlab.org\" target=\"_blank\"\u003ewww.vahdatlab.org\u003c/a\u003e\u003c/span\u003e\u003cspan address=\"http://www.vahdatlab.org\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. The main steps of the spinal cord processing pipeline are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e, including preprocessing, registration to the spinal cord template (PAM50 from SCT toolbox), and regression analysis, as described below.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003ePreprocessing\u003c/h2\u003e \u003cp\u003eFollowing removing of the first few volumes to reach a steady state condition (in our dataset, no volume was removed due to the inclusion of dummy scans during the acquisition), motion correction was performed on the EPI data in two steps (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). First, a 3d linear rigid-body transformation was performed using mcflirt tool (FSL) on the entire brain and spinal cord FOV. Next, the EPI was divided into the spinal cord and brain FOVs. Slice-wise (2d) non-linear motion correction was then performed in a relatively large (23 voxels diameter) cylindrical mask of the spinal cord using the \u003cem\u003esct_create_mask\u003c/em\u003e (SCT). The performance of various motion correction methods (3d linear alone, slice-wise nonlinear alone, or both) was calculated and compared using two indices: i) DVARS that measures changes in BOLD signal intensity from one volume to the next (Eq.\u0026nbsp;1, (Afyouni and Nichols, \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2018\u003c/span\u003e), and ii) tSNR (temporal signal-to-noise-ratio) which is the ratio of the mean by the standard deviation of the BOLD signal timeseries (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({Y}_{t}\\)\u003c/span\u003e\u003c/span\u003e) in each voxel (Eq.\u0026nbsp;2), \u003cem\u003en\u003c/em\u003e being the number of acquired volumes.\u003c/p\u003e \u003cp\u003eEquation 1\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\text{D}\\text{V}\\text{A}\\text{R}\\text{S}=\\sqrt{\\frac{1}{n}\\sum _{t=1}^{n}{({Y}_{t}-{Y}_{t-1})}^{2}}\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eEquation 2\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\text{t}\\text{S}\\text{N}\\text{R}= mean\\left({Y}_{t}\\right)/std\\left({Y}_{t}\\right)\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eFollowing motion correction, gradient-echo field map correction was performed using FUGUE tool (FSL) to correct for geometrical distortion due to magnetic field inhomogeneities (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). These artifacts are usually larger in the spinal cord compared to the brain due to the proximity to the thorax, vertebrae, and disks, and can cause a shift in the phase-encoding direction, which can be compensated at this step. Next, spatial smoothing was performed inside a mask of the spinal cord in a slice-wise fashion using \u003cem\u003e3dBlurToFWHM\u003c/em\u003e command from AFNI (Cox, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e1996\u003c/span\u003e). A kernel size of 3.2 mm (2 times the voxel size) was used for smoothing, although a 2.4 mm kernel size resulted in very similar task-based activation maps.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eRegistration to template\u003c/h2\u003e \u003cp\u003eThe spinal cord registration pipeline is displayed in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. Prior to registration, a spinal cord mask was generated for each of the T2w and T1w images using a deep learning segmentation algorithm as implemented in \u003cem\u003esct_deepseg_sc\u003c/em\u003e (SCT) with automatic centerline detection using the support vector machine and convolutional neural network options (Gros et al., \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Following visual inspection, manual corrections were carried out to address any minor inaccuracies in the spinal cord mask. Next, a centerline was drawn manually on the T2w and EPI images. To be consistent across participants, in the T2w image, the centerline in the dorsoventral direction was defined at a voxel location just posterior to the gray commissure, which was clearly visible in the acquired T2w image. For the EPI image, since the gray/white matter distinction is less clear, we used the calculated tSNR map to define the centerline. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, the cord exhibits noticeably higher tSNR values compared to the surroundings, including the cerebrospinal fluid (CSF), blood vessels, and connective tissues. Hence, we identified the centerline at the midpoint of the area with the highest tSNR (yellow area in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e; when tSNR map is thresholded at 40, color-coded in red-yellow).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eRegistration to the template (PAM50) was performed in three steps (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e): 1) EPI image registration to the T2w space using slice-wise centerline images alignment, 2) T2w image registration to the T1w space through the cord segmentation images using \u003cem\u003esct_register_multimodal\u003c/em\u003e (SCT), and 3) T1w image registration to PAM50 through cord segmentation using \u003cem\u003esct_register_to_template\u003c/em\u003e (SCT). In detail, for step 1, EPI and T2w images were registered by performing slice-wise 2d in-plane translations (in the dorsoventral and mediolateral directions) to align their centerlines using flirt command (FSL) with the nearest-neighbor interpolation. For step 2, we manually identified the C1-T1 disc labels in both T1w and T2w images, and used their calculated cord segmentations images to compute a slice-wise rigid body transformation using \u003cem\u003esct_register_multimodal\u003c/em\u003e command (any mismatch in the rostrocaudal direction was first corrected by aligning T1w and T2w disk label locations). For step 3, the cord segmentation and disc labels in the T1w image were employed for registration to PAM50 through \u003cem\u003esct_register_to_template\u003c/em\u003e command, using the default options (including straightening, vertebral alignment, and non-linear slice-wise registration). Note that only step 1 was applied to the 4d EPI images. The steps 2 and 3 warping transformations were concatenated using \u003cem\u003esct_concat_transfo\u003c/em\u003e, and were later applied to the summary statistics images (regression coefficients and their variance images) using \u003cem\u003esct_apply_transfo command\u003c/em\u003e.\u003c/p\u003e \u003cp\u003eSubsequently, we compared the performance of the proposed registration method with the gold-standard spinal cord registration approach which only utilizes the T1W image (Landelle et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). This approach was implemented by calculating the transformations from the EPI to T1W image using \u003cem\u003esct_register_multimodal\u003c/em\u003e, and from T1W image to PAM50 as calculated above. The two transformations are then concatenated to register the EPI image directly to PAM50. To evaluate the registration performance, the binary cord segmentation mask in the EPI space was registered to PAM50 space using each method, and the result was compared to the binary template cord mask. Two indices were defined by calculating: i) the number of voxels in the template mask that are missed in the registered mask, normalized by the total number of template cord voxels, and ii) the overlap between the registered and the template cord masks, normalized by the total number of voxels in both masks.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eRegression analysis\u003c/h2\u003e \u003cp\u003eWe performed participant-level general linear model (GLM) analyses on the preprocessed EPI data in the T2w space in two steps. First, a GLM was constructed to model and remove the effects of three sets of confounds regressors, including: i) the mean CSF signal averaged in a surround mask generated by dilating the mask of spinal cord (\u003cem\u003ei.e.\u003c/em\u003e, GM\u0026thinsp;+\u0026thinsp;WM), ii) 6 motion correction parameters (x, y, and z translations and rotations) derived from the linear motion correction step, and 2 slice-wise motion correction regressors generated by the non-linear motion correction step on the spinal cord FOV, and iii) 5 regressors, identified as confound, extracted from independent component analysis (ICA). The fastica algorithm (Hyvarinen, \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e1999\u003c/span\u003e) was used to decompose the preprocessed spinal cord EPI data (before spatial smoothing) within a large mask (including both the spinal cord and the surround mask) into 30 components (Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). We calculated 3 indicators to systematically sort components based on their likelihood of being task-related, as opposed to those related to physiological/scanner noise: 1) spatial indicator: the ratio of significantly active (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) voxels within the spinal cord mask to that within the surround mask, 2) frequency indicator: the ratio of the power of the component\u0026rsquo;s time-series frequency response around the task central frequency (in our block-design paradigm: 1 / (task period\u0026thinsp;+\u0026thinsp;rest period)\u0026thinsp;=\u0026thinsp;1/40 s\u0026thinsp;=\u0026thinsp;0.025 Hz; 0.01\u0026thinsp;\u0026lt;\u0026thinsp;\u003cem\u003ef\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.1 Hz may be used for resting-state data) to that outside the task-related frequency range with some margins (\u003cem\u003ef\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.025/1.5 and \u003cem\u003ef\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.025*1.5 Hz), and 3) temporal indicator: the Pearson correlation between the task timing signal (boxcar function in our task) convolved with a canonical hemodynamic response function (HRF) and the time-series of each component. A \u003cem\u003econfound index\u003c/em\u003e was defined by multiplying the 3 indicators; the lower the \u003cem\u003econfound index\u003c/em\u003e, the higher the likelihood that component represents a confound. The 5 components with the least \u003cem\u003econfound index\u003c/em\u003e were selected, and their time series were included in the first GLM as confound.\u003c/p\u003e \u003cp\u003eNext the residual (clean) image from the first GLM was high pass filtered at 1/60 Hz. The timing of the grip force block design task was modeled as a regressor of interest, while motion outliers were modeled as confounds. The latter were calculated using \u003cem\u003efsl_motion_outliers\u003c/em\u003e command on the non-preprocessed spinal cord EPI data (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The summary statistics images (\u003cem\u003eCope\u003c/em\u003e and \u003cem\u003eVarCope\u003c/em\u003e images in FSL) related to the regressor of interest were then registered to PAM50 space using the concatenated warping transformation calculated in the registration step.\u003c/p\u003e \u003cp\u003eWe generated a set of individual-specific masks to specify low tSNR voxels for each participant as follows. For each EPI run, a tSNR \u003cem\u003ethreshold\u003c/em\u003e was calculated using Maximum a Posteriori (MAP) estimation of two Gaussian distributions on the calculated tSNR histogram data inside a large spinal mask (including both the cord and the surround binary masks; Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). A Gaussian mixture distribution model using two components was fitted to the tSNR histogram using MATLAB \u003cem\u003efitgmdist\u003c/em\u003e function. The \u003cem\u003ethreshold\u003c/em\u003e was defined as the MAP estimate of the optimal boundary, which is the intersection of the two Gaussian distributions (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). Low tSNR voxels (tSNR\u0026thinsp;\u0026lt;\u0026thinsp;\u003cem\u003ethreshold\u003c/em\u003e) were identified for every EPI run and participant. A subject-specific low tSNR mask was then defined as the intersection of low tSNR voxels from all EPI runs for each participant (a total of 2 runs in our data). Note that the two Gaussian distributions represent, hypothetically, the spinal cord and the CSF/surround voxels, with high tSNR voxels expected to be located inside the spinal cord mask as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e. Hence, the goal of this step is to calculate an optimal threshold for differentiating spinal cord and CSF/surround voxels and identify those voxels inside the spinal cord mask that, inconsistently, possess low tSNR due to severe distortion and signal loss in the EPI spinal data.\u003c/p\u003e \u003cp\u003eFinally, a group-level GLM model was constructed to average individual-level maps using a mixed-effect (FLAME 1\u0026thinsp;+\u0026thinsp;2) model as implemented in FSL. The individual-specific low tSNR masks modeling the low tSNR voxels in each participant were entered as confound (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). First we generated a group-level spinal cord mask that contains voxels, in which the majority of participants showed an acceptable tSNR (\u003cem\u003ei.e.\u003c/em\u003e, tSNR\u0026thinsp;\u0026gt;\u0026thinsp;\u003cem\u003ethreshold\u003c/em\u003e). This condition ensures that for every voxel in the group-level mask, high quality data from at least half of the participants were included for averaging, and, at the same time, limits the number of confound regressors to half the degrees of freedom in the group-level GLM. Hence, a total of up to \u003cem\u003en\u003c/em\u003e = (number of participants / 2) 4d voxel-wise regressors were included to account for and remove the effect of low tSNR voxels in every participants from group-level averaging. This procedure enabled us to only incorporate voxels with an acceptable tSNR in the group-level averaging, without the need to remove these voxels from the group-level mask/analysis. All group-level statistical maps were then corrected for multiple comparisons using Gaussian random field (GRF) theory, \u003cem\u003eZ\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;2.7, and cluster-level threshold \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003eBrain fMRI processing\u003c/h2\u003e \u003cp\u003eBrain image preprocessing was performed with the FSL software package, using the same pipeline as described previously (Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). One participant was removed from brain analysis due to partial coverage of brain FOV. Preprocessing included the following steps: skull stripping (FSL, BET), motion correction (FSL, MCFLIRT), high-pass temporal filtering (1/60 Hz), GRE field map correction, spatial smoothing using an isotropic Gaussian kernel of 3.2 mm full width at half maximum. BBR method (FSL) was used to register the EPI data to the T1w anatomical image, followed by a nonlinear registration (FNIRT, FSL) from the T1w image to the MNI template. These two transformations were concatenated to register the summary statistics images from the participant level analysis in the EPI space to the MNI template. Participant-level GLM was constructed to identify brain areas activated during motor control task using a boxcar function modeling our block-design paradigm. The 6 motion correction parameters (x, y, and z translations and rotations) derived from the linear motion correction step were included as confound in this analysis. Furthermore, volumes detected as motion outliers (described earlier) were included as confound in the GLM. Next, a group-level GLM model was constructed to average participant-level maps using a mixed-effect (FLAME 1\u0026thinsp;+\u0026thinsp;2) model as implemented in FSL. All group-level statistical maps were then corrected for multiple comparisons using Gaussian random field (GRF) theory, \u003cem\u003eZ\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;2.7 (voxel-wise uncorrected threshold), and a cluster-wise corrected threshold at \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eFunctional connectivity analysis\u003c/h2\u003e \u003cp\u003eWe performed region of interest (ROI)-based functional connectivity analysis to examine the communication between the spinal cord and brain BOLD time-series. We defined 11 \u003cem\u003ea-priori\u003c/em\u003e selected ROIs in the brain shown to be part of the sensorimotor network using the MNI brain atlas and brainstem atlas (Cauzzo et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Singh et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). The brain ROIs included: left precentral gyrus, left postcentral gyrus, supplementary motor area, left thalamus, right cerebellum lobule I-IV, right cerebellum lobule V-VI, left red nucleus, bilateral pontine reticular formation, bilateral medullary reticular formation. We also defined 2 ROIs in the spinal cord (right ventral horn C5-C8, right dorsal horn C5-C8 segmental levels). To restrict our analysis to subregions somatotopically related to the forelimb region, each ROI only included the most active voxels during handgrip force production by selecting voxels with the highest \u003cem\u003eZ\u003c/em\u003e-score from the individual-level GLM analysis (top 10, 20, 30, or 40th percentile voxels within each ROI). We then calculated the average clean BOLD signal (the residual image from the first GLM) within each ROI at each percentile threshold. For each participant and run, the partial correlation (Marrelec et al., \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) was calculated between each pair of the brain and spinal cord ROI\u0026rsquo;s time-series (22 pairs), by accounting for and removing the effects of task presentation (task timing convolved with HRF). The average partial correlation across 4 percentile thresholds were calculated and reported as the functional connectivity (FC) for each ROI pair. For each participant, FC values between the two runs were averaged and Fisher transformation was applied to meet normal distribution assumption. Finally, \u003cem\u003et\u003c/em\u003e-statistics were performed to identify significant functional connections during the motor task between the brain and spial cord ROIs across participants (corrected for multiple comparisons using Bonferroni correction, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eResting-state analysis\u003c/h2\u003e \u003cp\u003e \u003cb\u003eDenoising pipelines.\u003c/b\u003e To investigate the impact of different sources of noise we compared the performance of different denoising methods. For each methods, the fMRI time-series was regressed with various variables of no-interest (confounds) using FSL's fMRI Expert Analysis Tool (FEAT). No filtering was applied at this step. All methods incorporated 8 motion-related regressors, 6 obtained from the linear motion correction step, and 2 slice-wise regressors from the SCT motion correction step, as described earlier. Next, we extracted the mean signal of the CSF in a surround mask generated by dilating the mask of the spinal cord. For each participant, the physiological noise modelling was calculated using the Tapas Physio toolbox (Kasper et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). We used the RETROspective Image CORrection (RETROICOR) procedure (Glover et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2000\u003c/span\u003e) to generate noise regressors from peripheral physiological recordings (heart rate and respiration). Specifically, we modeled, four cardiac, four respiratory, harmonics, and four multiplicative terms for the interactions between respiratory and cardiac noise (32 regressors in total, similar to (Brooks et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Kaptan et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Kong et al., \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). To identify non-neural fluctuations, we extracted the first five principal components of the unfiltered and unsmoothed CSF signal (in the surround mask) in the participant\u0026rsquo;s native space using the CompCor approach (Behzadi et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2007\u003c/span\u003e). Additional regressors, identified as confound, were extracted from ICA in a dilated mask of the spinal cord as described in the \u003cspan refid=\"Sec9\" class=\"InternalRef\"\u003eRegression analysis\u003c/span\u003e section (3.89\u0026thinsp;\u0026plusmn;\u0026thinsp;1.66 components).\u003c/p\u003e \u003cp\u003eFinally, a specific set of regressors of no-interest was accounted for in separate GLMs, corresponding to different denoising models as follow:\u003c/p\u003e \u003cp\u003e \u003col\u003e \u003cspan\u003e \u003cli\u003e \u003cp\u003e\u0026ldquo;\u003cem\u003eMeanCSF\u003c/em\u003e\u0026rdquo; (baseline parameters, mean CSF)\u003c/p\u003e \u003c/li\u003e \u003c/span\u003e \u003cspan\u003e \u003cli\u003e \u003cp\u003e\u0026ldquo;\u003cem\u003eICA\u003c/em\u003e\u0026rdquo; (baseline parameters, mean CSF, ICA components)\u003c/p\u003e \u003c/li\u003e \u003c/span\u003e \u003cspan\u003e \u003cli\u003e \u003cp\u003e\u0026ldquo;\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; (baseline parameters, mean CSF, 5 CompCor components)\u003c/p\u003e \u003c/li\u003e \u003c/span\u003e \u003cspan\u003e \u003cli\u003e \u003cp\u003e\u0026ldquo;\u003cem\u003eRetroicor\u003c/em\u003e\u0026rdquo; (baseline parameters, mean CSF, 32 Retroicor components)\u003c/p\u003e \u003c/li\u003e \u003c/span\u003e \u003cspan\u003e \u003cli\u003e \u003cp\u003e\u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; (baseline parameters, mean CSF, ICA and 5 CompCor components)\u003c/p\u003e \u003c/li\u003e \u003c/span\u003e \u003c/ol\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eVariance of the model.\u003c/b\u003e For each method, we calculated the voxel-wise explained variance of the model (R\u003csup\u003e2\u003c/sup\u003e) as a performance metric to provide deeper insights into the impact of noise removal from various sources using the following equation:\u003c/p\u003e \u003cp\u003eEquation 3\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{R}}^{2}=[var\\left({Y}_{t}\\right) - var\\left({\\epsilon }_{t}\\right)]/var\\left({Y}_{t}\\right)\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eWhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(var\\left({Y}_{t}\\right)\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(var\\left({\\epsilon }_{t}\\right)\\)\u003c/span\u003e\u003c/span\u003e denote respectively the variances of the fMRI times-series and the residual signal (\u003cem\u003eres4d.nii\u003c/em\u003e in FSL) at each voxel, after applying the regression model associated with each method. Next, the average R\u003csup\u003e2\u003c/sup\u003e was calculated inside the surround mask (\u003cem\u003ei.e.\u003c/em\u003e, the CSF mask), as well as the spinal cord mask. The ratio between the R\u003csup\u003e2\u003c/sup\u003e of CSF and spinal cord (R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003eCSF\u003c/sub\u003e / R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003ecord\u003c/sub\u003e) is reported to evaluate the model's performance. A good model for removing the physiological noise would account for more variance in the CSF, while for less variance related to neural activation in the spinal cord. Hence, the larger the R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003eCSF\u003c/sub\u003e / R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003ecord\u003c/sub\u003e ratio, the greater the capacity of the model in explaining the BOLD signal variability inside the CSF compared to the spinal cord, and the better the model in removing the physiological noise.\u003c/p\u003e \u003cp\u003eBeing aware of the limitations of the R\u003csup\u003e2\u003c/sup\u003e as a quality check metric, such as its failure to account for each model's degree of freedom, we also calculated the adjusted variance of the model (AdjR\u003csup\u003e2\u003c/sup\u003e) as follows:\u003c/p\u003e \u003cp\u003eEquation 4\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{A}\\text{d}\\text{j}\\text{R}}^{2}=1- \\frac{(1-{\\text{R}}^{2})(N-1)}{N-P-1}\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eWhere \u003cem\u003eN\u003c/em\u003e is the number of fMRI image time points, and \u003cem\u003eP\u003c/em\u003e corresponds to the number of regressors of no-interests (confounds) included in each model. Unlike the R\u003csup\u003e2\u003c/sup\u003e, which measures the proportion of variance explained by the model, AdjR\u003csup\u003e2\u003c/sup\u003e takes into account the degree of freedom of the model. This implies that with an increase in the number of regressors (resulting in the degree of freedom increase), the AdjR\u003csup\u003e2\u003c/sup\u003e value will increase only if the additional regressors enhance the overall model performance. The AdjR\u003csup\u003e2\u003c/sup\u003e was extracted in the CSF mask. t-statistics were performed to compare the R\u003csup\u003e2\u003c/sup\u003e and AjdR\u003csup\u003e2\u003c/sup\u003e between each model, and differences were considered significant with a \u003cem\u003ep\u003c/em\u003e-value\u0026thinsp;\u0026lt;\u0026thinsp;0.005, which was adjusted for multiple comparisons using Bonferroni correction (0.05 /10\u0026thinsp;=\u0026thinsp;0.005).\u003c/p\u003e \u003cp\u003e \u003cb\u003eResting-state functional connectivity.\u003c/b\u003e To assess the impact of denoising methods on functional connectivity measures, we conducted ROI-based analyses on the residual (cleaned) images. We defined 4 ROIs in the spinal cord (right ventral horn C5-C8, right dorsal horn C5-C8, left ventral horn C5-C8, left dorsal horn C5-C8) and 1 ROI in the surround mask cord (CSF C5-C8). We extracted and band-pass filtered (0.01\u0026ndash;0.17 Hz)(Barry et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2018a\u003c/span\u003e) the time series within each ROI and calculated the average signal. Next, for each participant, we computed the partial correlation between each pair of spinal cord ROIs while the CSF average signal was used as a covariate. The Fisher-transformed correlation coefficients were employed to establish functional connectivity between each ROI pair, specifically between ventral horns (Ventro-Ventral FC), dorsal horns (Dorso-Dorsal DC), ventro-dorsal right horns (Right FC) and ventro-dorsal left horns (Left FC). For each FC, t-statistics were performed to compare the functional connectivity between each model (corrected for multiple comparisons using Bonferroni correction 0.05 /10\u0026thinsp;=\u0026thinsp;0.005).\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003eWe collected simultaneous spinal cord-brain fMRI data, covering the whole brain and the entire cervical cord down to the first thoracic vertebra, using simultaneous multi-slice (SMS) EPI sequence, together with the application of neck and brachial plexus SatPad\u0026trade; set, as described in the \u003cspan refid=\"Sec2\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section. This image acquisition protocol provided a reasonable spatial and temporal resolution (1.6x1.6 mm\u003csup\u003e2\u003c/sup\u003e in plane, TR: 1.85 ms), as well as a reasonable tSNR in the spinal cord (22.32\u0026thinsp;\u0026plusmn;\u0026thinsp;0.86 mean\u0026thinsp;\u0026plusmn;\u0026thinsp;sem; Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e) of our group of participants. Furthermore, the low tSNR voxels were statistically excluded, per participant, from the group-level analysis. On average 76\u0026thinsp;\u0026plusmn;\u0026thinsp;2.1 (mean\u0026thinsp;\u0026plusmn;\u0026thinsp;sem) percent of the spinal cord voxels passed the MAP threshold and were included in group-level analysis. The average tSNR in the voxels included in the group-level analysis was 25.40\u0026thinsp;\u0026plusmn;\u0026thinsp;0.81 (mean\u0026thinsp;\u0026plusmn;\u0026thinsp;sem) across participants (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e).\u003c/p\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eMotion correction procedures\u003c/h2\u003e \u003cp\u003eFirst, we assessed the efficiency of various motion correction procedures in the spinal cord using DVARS (Eq.\u0026nbsp;1) and the spinal cord tSNR metrics (Eq.\u0026nbsp;2). Optimal outcomes encompass reduced DVARS and enhanced spinal cord tSNR values. The performance of 3 methods were compared to the raw (no motion correction) data, including 3d linear motion correction on the entire brain and spinal cord FOV (SP\u0026thinsp;+\u0026thinsp;BR Flirt), slice-wise nonlinear correction on the spinal cord FOV (Slice-wise SCT), and the combined procedure (first performing 3d linear correction on the entire FOV, followed by slice-wise nonlinear correction on the spinal cord slices) as implemented in FASB. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e, the combined procedure outperforms the other methods in terms of DVARS and tSNR (repeated-measures ANOVA, corrected for multiple comparisons using Bonferroni correction, \u003cb\u003eDVARS\u003c/b\u003e: \u003cem\u003eF(3,42)\u003c/em\u003e\u0026thinsp;=\u0026thinsp;14.3, \u003cem\u003ep\u003c/em\u003e (FASB vs Raw)\u0026thinsp;\u0026lt;\u0026thinsp;0.005, \u003cem\u003ep\u003c/em\u003e (FASB vs Flirt)\u0026thinsp;\u0026lt;\u0026thinsp;0.05, \u003cem\u003ep\u003c/em\u003e (FASB vs SCT)\u0026thinsp;\u0026lt;\u0026thinsp;0.05, \u003cem\u003ep\u003c/em\u003e (SCT vs Flirt)\u0026thinsp;\u0026gt;\u0026thinsp;0.05; \u003cb\u003eSpinal cord tSNR\u003c/b\u003e: \u003cem\u003eF(3,42)\u003c/em\u003e\u0026thinsp;=\u0026thinsp;23.5, \u003cem\u003ep\u003c/em\u003e (FASB vs Raw)\u0026thinsp;\u0026lt;\u0026thinsp;0.001, \u003cem\u003ep\u003c/em\u003e (FASB vs Flirt)\u0026thinsp;\u0026lt;\u0026thinsp;0.01, \u003cem\u003ep\u003c/em\u003e (FASB vs SCT)\u0026thinsp;\u0026lt;\u0026thinsp;0.01, \u003cem\u003ep\u003c/em\u003e (SCT vs Flirt)\u0026thinsp;\u0026gt;\u0026thinsp;0.05). This analysis revealed that combining both 3d linear and slice-wise non-linear motion correction procedures yields reduced temporal variability and enhanced temporal SNR in the spinal cord.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eRegistration procedures\u003c/h2\u003e \u003cp\u003eNext, we compared the performance of various registration methods to align the EPI spinal cord image to a template (PAM50). As outcome measures, we calculated percentage error and overlap between the registered binary EPI spinal cord mask and the PAM50 spinal cord segmentation, as reported in Methods. Smaller \u003cem\u003eError\u003c/em\u003e and larger \u003cem\u003eOverlap\u003c/em\u003e percentage values are desirable. As described in the Methods, the performances of 3 procedures are compared, including registration using T1W image as the intermediate image (T1W-based), GRE field map correction together with T1W-based registration (GRE\u0026thinsp;+\u0026thinsp;T1W-based), and FASB registration procedure as summarized in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, the FASB registration procedure outperforms the other registration methods (which do not include T2w image) in terms of both \u003cem\u003eError\u003c/em\u003e and \u003cem\u003eOverlap\u003c/em\u003e for the spinal cord masks (repeated-measures ANOVA, corrected for multiple comparisons using Bonferroni correction, \u003cb\u003eError\u003c/b\u003e: \u003cem\u003eF(2,28)\u003c/em\u003e\u0026thinsp;=\u0026thinsp;37, \u003cem\u003ep\u003c/em\u003e (FASB vs T1W)\u0026thinsp;\u0026lt;\u0026thinsp;0.001, \u003cem\u003ep\u003c/em\u003e (FASB vs GRE\u0026thinsp;+\u0026thinsp;T1W)\u0026thinsp;\u0026lt;\u0026thinsp;0.001, \u003cem\u003ep\u003c/em\u003e (T1W vs GRE\u0026thinsp;+\u0026thinsp;T1W)\u0026thinsp;\u0026gt;\u0026thinsp;0.05; \u003cb\u003eOverlap\u003c/b\u003e: \u003cem\u003eF(2,28)\u003c/em\u003e\u0026thinsp;=\u0026thinsp;11.5, \u003cem\u003ep\u003c/em\u003e (FASB vs T1W)\u0026thinsp;\u0026lt;\u0026thinsp;0.001, \u003cem\u003ep\u003c/em\u003e (FASB vs GRE\u0026thinsp;+\u0026thinsp;T1W)\u0026thinsp;\u0026lt;\u0026thinsp;0.05, \u003cem\u003ep\u003c/em\u003e (T1W vs GRE\u0026thinsp;+\u0026thinsp;T1W)\u0026thinsp;\u0026gt;\u0026thinsp;0.05). This analysis shows that performing the additional slice-wise centerline alignment between EPI and T2w images, as proposed in FASB, significantly improves the overall accuracy of the spinal cord registration to the template space (from 78\u0026ndash;84% overlap).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003ePhysiological noise reduction procedures\u003c/h2\u003e \u003cp\u003eNext, we compared the performance of different physiological noise removal models on the collected resting-state fMRI data. We compared 5 models (see Methods for details): \u003cem\u003e\u0026ldquo;meanCSF\u0026rdquo;\u003c/em\u003e (the baseline model for removing motion correction parameters and the CSF average), \u003cem\u003e\u0026ldquo;ICA\u0026rdquo;\u003c/em\u003e (as implemented in FASB), \u003cem\u003e\u0026ldquo;CompCor\u0026rdquo;\u003c/em\u003e (a PCA-based model), \u003cem\u003e\u0026ldquo;Retroico\u003c/em\u003er\u003cem\u003e\u0026rdquo;\u003c/em\u003e (based on the recorded cardiac and respiratory signals), and \u003cem\u003e\u0026ldquo;ICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo;\u003c/em\u003e (combined \u003cem\u003e\u0026ldquo;ICA\u0026rdquo;\u003c/em\u003e and \u003cem\u003e\u0026ldquo;CompCor\u0026rdquo;\u003c/em\u003e models). First, we compared the ratio between the R\u003csup\u003e2\u003c/sup\u003e of CSF and spinal cord (R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003eCSF\u003c/sub\u003e / R\u003csup\u003e2\u003c/sup\u003e\u003csub\u003ecord\u003c/sub\u003e) between all denoising models to evaluate model performances (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). A higher value indicates that the variance of the CSF is explained more substantially by the model relative the variance of the cord. Results showed a significant increase across most of the denoising models compared to the baseline model \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo; (\"\u003cem\u003eICA\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e: t(17)\u0026thinsp;=\u0026thinsp;4.67, p\u0026thinsp;=\u0026thinsp;0.0002; \"\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u0026rdquo;\u003c/em\u003e: t(17)\u0026thinsp;=\u0026thinsp;4.87, p\u0026thinsp;=\u0026thinsp;0.0001; \"\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e: t(17)\u0026thinsp;=\u0026thinsp;4.34, p\u0026thinsp;=\u0026thinsp;0.0004), with the exception of the \u0026ldquo;\u003cem\u003eRetroico\u003c/em\u003er\u0026rdquo; model where no significant difference was observed with the \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo; model (t(17)\u0026thinsp;=\u0026thinsp;\u003cem\u003e1.12, p\u0026thinsp;\u0026lt;\u0026thinsp;0.27\u003c/em\u003e). In addition, \u0026ldquo;\u003cem\u003eRetroico\u003c/em\u003er\u0026rdquo; model showed a lower ratio than \u0026ldquo;CompCor\u0026rdquo; (t(17)= -5.38, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001), \u0026ldquo;ICA\u0026rdquo; (t(17)= -4.02, p\u0026thinsp;=\u0026thinsp;0.0009) and \u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; (t(17)\u0026thinsp;=\u0026thinsp;5.97, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001) models. \u003cem\u003e\u0026ldquo;CompCor\u0026rdquo;\u003c/em\u003e model had significant higher ratio than \u003cem\u003e\u0026ldquo;ICA\u0026rdquo; (\u003c/em\u003e\"\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003eICA\u0026rdquo; (t(17)\u0026thinsp;=\u0026thinsp;3.45, p\u0026thinsp;=\u0026thinsp;0.003)\u003c/em\u003e, while no significant differences, after the Bonferroni correction, were observed between the \u003cem\u003e\u0026ldquo;ICA\u003c/em\u003e\u0026rdquo; and \u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo;\u003c/em\u003e (t(17)\u0026thinsp;=\u0026thinsp;3.03, p\u0026thinsp;=\u0026thinsp;0.007\u003cem\u003e)\u003c/em\u003e, or \"\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; and \u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo; (t(17)\u0026thinsp;=\u0026thinsp;2.60, p\u0026thinsp;=\u0026thinsp;0.018\u003c/em\u003e).\u003c/p\u003e \u003cp\u003eTo account for the degrees of freedom of each model (depending on the number of confound regressors), we also calculated the adjusted R\u003csup\u003e2\u003c/sup\u003e (AdjR\u003csup\u003e2\u003c/sup\u003e; Eq.\u0026nbsp;4) inside the CSF mask. Adding new confound regressors does increase the AdjR\u003csup\u003e2\u003c/sup\u003e value only if the additional regressors enhance the overall model performance. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB, all models have higher AdjR\u003csup\u003e2\u003c/sup\u003e than the baseline model (\u003cem\u003e\u0026ldquo;ICA\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo;: t(17)\u0026thinsp;=\u0026thinsp;12.89, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; \u003cem\u003e\u0026ldquo;CompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo; : t(17)\u0026thinsp;=\u0026thinsp;24.25, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; \u003cem\u003e\u0026ldquo;Retroicor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u0026rdquo;: t(17)\u0026thinsp;=\u0026thinsp;7.01, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; \u0026ldquo;ICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo;: t(17)\u0026thinsp;=\u0026thinsp;24.38, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001). \u0026ldquo;\u003cem\u003eRetroico\u003c/em\u003er\u0026rdquo; model showed lower AdjR\u003csup\u003e2\u003c/sup\u003e than \u0026ldquo;\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; (t(17)= -20.40, p\u0026thinsp;=\u0026thinsp;0.0001), \u0026ldquo;ICA\u0026rdquo; (t(17)= -7.07, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001) and \u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo;\u003c/em\u003e (t(17)= -19.45, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001) models. \u003cem\u003e\u0026ldquo;ICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo;\u003c/em\u003e model had significantly higher AdjR\u003csup\u003e2\u003c/sup\u003e than \u003cem\u003e\u0026ldquo;ICA\u0026rdquo;\u003c/em\u003e (t(17)\u0026thinsp;=\u0026thinsp;20.40, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001) and \u003cem\u003e\u0026ldquo;CompCor\u0026rdquo;\u003c/em\u003e models (t(17)\u0026thinsp;=\u0026thinsp;9.33, p\u0026thinsp;\u0026lt;\u0026thinsp;0001). Finally \u0026ldquo;\u003cem\u003eCompCor\u0026rdquo;\u003c/em\u003e model showed significantly higher AdjR\u003csup\u003e2\u003c/sup\u003e than \u003cem\u003e\u0026ldquo;ICA\u0026rdquo;\u003c/em\u003e model (t(17)\u0026thinsp;=\u0026thinsp;10.35, p\u0026thinsp;\u0026lt;\u0026thinsp;0001).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eAssessing the impact of denoising procedures on functional connectivity\u003c/h2\u003e \u003cp\u003eNext to assess the impact of denoising methods on functional connectivity (FC) measures, we conducted ROI-based analyses (see Methods). As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC, no significant differences in FC between the different denoising models were observed associated with ipsilateral functional connectivity measures (ventro-dorsal right, and ventro-dorsal left). \u003cem\u003e\u0026ldquo;ICA\u0026thinsp;+\u0026thinsp;CompCor\u0026rdquo;\u003c/em\u003e model led to significantly reduced functional connectivity between the left and right ventral horns compared to most of the other models (\u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo;: t(17)= -3.86, p\u0026thinsp;=\u0026thinsp;0.001; \u0026ldquo;\u003cem\u003eICA\u003c/em\u003e\u0026rdquo;: t(17)= -3.80; p\u0026thinsp;=\u0026thinsp;0.001; \u0026ldquo;\u003cem\u003eRetroicor\u003c/em\u003e\u0026rdquo;: t(17)= -4.36, p\u0026thinsp;=\u0026thinsp;0004) and a trend toward significant when compared to the \u0026ldquo;\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; model (\u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; \u003cem\u003evs\u003c/em\u003e. \u0026ldquo;\u003cem\u003eCompCor\u003c/em\u003e\u0026rdquo; t(17)\u0026thinsp;=\u0026thinsp;3.1, p\u0026thinsp;=\u0026thinsp;0.006). Similarly, the FC between the left and right dorsal horns was significantly lower for the \u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo; model compared to the others (\u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo;: t(17)= -5.59, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001, \u0026ldquo;\u003cem\u003eICA\u003c/em\u003e\u0026rdquo;: t(17)= -4.02, p\u0026thinsp;=\u0026thinsp;0.0008, \u0026ldquo;\u003cem\u003eCompcor\u003c/em\u003e\u0026rdquo;: t(17)= -3.62, p\u0026thinsp;=\u0026thinsp;0.002, \u0026ldquo;\u003cem\u003eRetroicor\u003c/em\u003e\u0026rdquo;: t(17)= -6.19, p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001). The FC was also found lower with \u0026ldquo;ICA\u0026rdquo; model compared to \u0026ldquo;\u003cem\u003emeanCSF\u003c/em\u003e\u0026rdquo; (t(17):-4.46, p\u0026thinsp;=\u0026thinsp;0.0003) and \u0026ldquo;\u003cem\u003eRetroicor\u003c/em\u003e\u0026rdquo; (t(17)= -4.56, p\u0026thinsp;=\u0026thinsp;0.0002) models. In sum, ventro-ventral and dorso-dorsal connections showed reduction in functional connectivity in the more complex model (\u0026ldquo;\u003cem\u003eICA\u0026thinsp;+\u0026thinsp;CompCor\u003c/em\u003e\u0026rdquo;) compared to the others. Interestingly for the weaker FC observed within horns (ventro-dorsal FC), the various denoising models did not significantly impact the FC.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eTask-related results\u003c/h2\u003e \u003cp\u003eFigure \u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e shows activation clusters associated with the motor task performance in the brain and brainstem. At the brain level, all of the anticipated sensorimotor cortical and subcortical areas related to unilateral hand movement (Vahdat et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2017\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2011\u003c/span\u003e) showed significant activation, including the contralateral primary motor cortex (M1), dorsal premotor cortex, primary and secondary somatosensory cortices (SI, SII), thalamus (ventral lateral and pulvinar nuclei), putamen, ipsilateral cerebellar cortex (lobules I-IV, and V-VI), bilateral posterior parietal cortex (PPC), and supplementary motor area (SMA) (corrected for multiple comparisons using GRF, \u003cem\u003eZ\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;2.7, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.5, Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). Interestingly, our high in-plane spatial resolution (1.6\u0026times;1.6 mm\u003csup\u003e2\u003c/sup\u003e) allowed us to detect and delineate task-related brainstem activations, including bilateral pontine reticular formation (PRF), medullary reticular formation (MRF), and red nucleus (RN) (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eB). Consistent with the laterality of descending projections from the reticular formation and RN (Davidson et al., \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Lemon, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Soteropoulos et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2012\u003c/span\u003e), task-related activation was greater in the ipsilateral MRF and PRF, and the contralateral RN. These results show that the acquired EPI brain images in our simultaneous protocol not only provide enough power to replicate previous task-related brain networks, but also offer detection and delineation of activity in important brainstem nuclei.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eNext, we also compared task-related activation maps in the spinal cord using different approaches. Four methods were implemented and compared as described in the \u003cspan refid=\"Sec2\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section: i) regular group-level regression, without the inclusion of ICA confounds (NoICA-NotSNR), ii) regular group-level regression, with the inclusion of ICA confounds (ICA-NotSNR), iii) group-level regression by inclusion of voxel-wise low-tSNR confounds, but without ICA confounds (NoICA-tSNR), and iv) group-level regression by inclusion of voxel-wise low-tSNR confounds, as well as ICA confounds (ICA-tSNR, implemented in FASB). As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e, all 4 methods resulted in significant activation clusters at the C5-C8 cervical segmental levels (corrected for multiple comparisons using GRF, \u003cem\u003eZ\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;2.7, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Compared to the other methods, FASB (ICA-tSNR) resulted in the largest number of activation clusters (4 clusters, including C5, C6, C7, and C8) and detected activation across all related cervical levels for the control hand/wrist muscles (Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Importantly, the largest peaks of activation were unequivocally located in the ipsilateral ventral horn, where the hand/wrist motoneurons are located. The activation at the C6-C7 level was largely restricted to the ventral horn, while it included the dorsal horn at the C8 segmental level. These results show that the acquisition/analysis pipeline proposed in FASB can differentiate task-related spinal cord activation based on myotome/dermatome (spinal level), side (left/right), and sensory-motor (dorsal/ventral) characteristics.\u003c/p\u003e \u003cp\u003eTo evaluate the connectivity related to the ascending and descending projections, we next investigated the partial correlation between \u003cem\u003ea-priori\u003c/em\u003e selected set of ROIs in the brain and spinal cord, when the effects of task paradigm are removed from this relationship (see Methods for details). The right (ipsilateral) ventral horn at C5-C8 cervical level showed significant functional connectivity to the left postcentral gyrus, SMA, thalamus, and right medullary reticular formation and cerebellum lobule V-VI (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e; \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Also, the right dorsal horn at C5-C8 cervical level showed significant functional connectivity to the left precentral and postcentral gyrus, SMA, thalamus, and right cerebellum lobule V-VI. These results demonstrate that significant task-related functional connectivity between the brain and spinal cord ROIs can be identified using FASB, consistent with the neuroanatomy of well-established descending and ascending projections in humans (Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Lemon, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2008\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eWe developed and validated an integrated processing pipeline for functional analysis of simultaneous spinal cord-brain EPI data (FASB), together with recommendations on collecting appropriate anatomical and functional data for that pipeline. The proposed EPI acquisition allowed us to acquire spinal cord images using a Siemens 3T MRI scanner, with low distortion due to magnetic susceptibility artifacts, and with high SNR (mean temporal SNR\u0026thinsp;\u0026gt;\u0026thinsp;20 in the cord), while having a relatively high spatial resolution at the spinal cord level (1.6 mm in-plane), and a reasonable acquisition time (TR\u0026thinsp;=\u0026thinsp;1.85 s). The novel aspects of the FASB pipeline are threefold (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e): 1) Improved motion correction of spinal cord fMRI volumes using both 3d linear and slice-wise non-linear transformations, 2) Improved spinal cord EPI registration to the template space (PAM50) using two intermediate anatomical images (T2w and T1w), and slice-wise centerline alignment to T2w image guided by the tSNR map of the EPI data (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), and 3) Improved detection power of the group-level analysis by removing the effects of low tSNR voxels in each participant, typically at the disk level, based on the MAP estimate of the tSNR threshold between the spinal cord and CSF/surround voxels (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eThe simultaneous acquisition of brain and spinal cord enabled us to use spatial information from both FOVs for motion correction of the spinal cord EPI volumes. Based on our results reported in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e, performing a linear transformation on the entire brain and spinal cord FOVs, followed by slice-wise non-linear correction at each spinal cord slice significantly improves the overall motion correction performance. This is likely due to the 3d linear correction step, which effectively addresses motion correction along the rostro-caudal direction, where the slice-wise correction step is insensitive.\u003c/p\u003e \u003cp\u003eFurthermore, FASB registration of the EPI to the template space significantly improves the spinal cord registration compared to current spinal registration approaches (Landelle et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Spinal cord EPI images are often spatially distorted at the disk level, and performing a nonlinear transformation generates non-optimal twisted warping fields. Magnetic susceptibility artifacts also result in signal loss in the distorted voxels, which makes the inclusion of these voxels in further analyses even more problematic. Hence, we proposed performing a simple slice-wise centerline alignment to match the centerlines of the EPI and T2w images. This is possible because we acquire a T2w image at the same angle/slice thickness as the EPI image. Moreover, to standardize the selection of centerline specifically in the distorted slices, we propose to select the EPI centerline based on the estimated tSNR map (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). Lastly, we account for and remove the effects of distorted slices later by including the signal of participant-specific low tSNR voxels as confound in the group-level regression analysis. Although the T2w medic image is one of the requirements for our registration pipeline, because of its sharp contrast for discerning the gray and white matter, it can also provide useful morphometric information that may be of interest in different clinical populations.\u003c/p\u003e \u003cp\u003eFor systematic selection of the tSNR threshold between the spinal cord and CSF/surround voxels, we propose to use the optimal boundary (MAP estimate) between two Gaussian distributions, fitted on the tSNR histogram. This makes our threshold selection procedure unbiased and dependent upon the quality of the acquired functional data. In our collected EPI data, the average tSNR threshold was 12.96 (\u0026plusmn;\u0026thinsp;0.40 sem). This allowed us to identify the \u0026ldquo;low tSNR voxels\u0026rdquo; (\u003cem\u003ei.e.\u003c/em\u003e, those inside the spinal cord mask whose tSNR was lower than the MAP estimate), and remove their effect at the group-level analysis by the inclusion of participant-specific voxel-wise regressors. To the best of our knowledge, our proposed method provides the only available solution to account for and remove the effects of low tSNR voxels from the group-level results, which is particularly critical in spinal cord fMRI due to large magnetic susceptibility artifacts at the spinal cord level.\u003c/p\u003e \u003cp\u003eFor removal of the scanner noise/physiological artifacts, which are particularly intensified in spinal cord fMRI, we propose to use ICA together with automatic selection of confound components based on their spatial, frequency, and temporal characteristics. Our results show that the performance of ICA and \u003cem\u003eRETROICOR\u003c/em\u003e (Glover et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2000\u003c/span\u003e) (based on recorded cardiac and respiratory rates) methods are comparable in terms of reducing the variability attributed to noise (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA-B). Combining the confound regressors extracted using \u003cem\u003eCompCor\u003c/em\u003e method (Behzadi et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2007\u003c/span\u003e) (a PCA-based approach) and ICA improves the noise removal performance as evaluated by the explained variance in CSF. However, it should be noted that CSF outcome measures are generally biased toward \u003cem\u003eCompCor\u003c/em\u003e, as this method performs PCA inside the CSF mask. Importantly, the analysis of functional connectivity between various spinal cord quadrants (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC) revealed that the various denoising models (\u003cem\u003eCompCor\u003c/em\u003e, ICA, \u003cem\u003eRETROICOR\u003c/em\u003e) did not significantly impact the calculated functional connectivity results. Hence, we included ICA as the method of choice in FASB, due to the ease of acquisition (no need for cardiac and respiratory recordings). However, note that our ICA confound selection procedure is particularly suited for task-based fMRI, in which the frequency and temporal criteria can be defined based on the task timing. In resting-state paradigms selection of confound components can still be achieved (e.g., more active voxels inside the CSF than the spinal cord mask, and larger power outside the neural-related frequency range [0.01 Hz, 0.17 Hz] (Barry et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2018a\u003c/span\u003e) than within this range), though it is not as straightforward as in the bask-based design. Hence, for resting-state analysis, a combination of ICA, \u003cem\u003eCompCor\u003c/em\u003e or \u003cem\u003eRETROICOR\u003c/em\u003e methods may be considered.\u003c/p\u003e \u003cp\u003eFor validation purposes, we examined the activation maps in the brain and spinal cord using FASB during the performance of a hand force control task. The force control task was selected due to 3 main advantages compared to similar motor tasks performed in the scanner: 1) The task is isometric, so there is no difference in movement kinematics across participants, 2) The produced force is normalized across individuals using their MVC, and 3) The isometric task minimizes motion artifacts, which are critical for spinal cord imaging at the neck level. The activated cortical and subcortical areas during this task are consistent with our previous fMRI studies using hand movements (Braa\u0026szlig; et al., \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2017\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2011\u003c/span\u003e), including contralateral M1, S1, premotor cortex, SII, PPC, thalamus, striatum, SMA, and ipsilateral cerebellum (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e). Importantly, we found significant task-related activations at the C5-C8 cervical level, mainly at the ipsilateral ventral horn (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e), consistent with the known neuroanatomical location of the finger/wrist motoneurons. Compared to our previous simultaneous spinal cord-brain fMRI studies (Khatibi et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Kinany et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2015\u003c/span\u003e), in which inclusion of more than 25 participants were typically recruited to detect significant spinal cord activations, analysis of only 15 participants using the FASB pipeline resulted in significant focal spinal cord activation.\u003c/p\u003e \u003cp\u003eImportantly, simultaneous spinal cord-brain fMRI enables investigation of functional connectivity related to the descending and ascending pathways (Landelle et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Our ROI-based analyses demonstrated significant functional connectivity during task performance between cortical, thalamic, cerebellar, and bulbar regions with the cervical ventral horn, and between cortical, thalamic and cerebellar regions with the cervical dorsal horn (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e). Interestingly, we found significant functional connectivity between the ipsilateral medullary reticular formation and the cervical ventral horn, consistent with the involvement of the reticulospinal pathway in hand motor control (Baker, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Soteropoulos et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). Since the goal of the current study was to validate FASB pipeline for spinal cord fMRI processing, we performed a simple ROI-based partial correlation analysis to evaluate spinal cord-brain functional connectivity. Partial correlation analysis can account for and remove the effects of task presentation from the BOLD signal time-series, so that their correlation does not simply represent co-activation of areas during task periods as compared to rest (Marrelec et al., \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2006\u003c/span\u003e). Other methods based on ICA (Vahdat et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2020\u003c/span\u003e, \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2012\u003c/span\u003e), or BASCO (Rissman et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2004\u003c/span\u003e) (Beta series correlation, particularly suited for an event-related task design) may be used to assess functional connectivity between the brain and spinal cord during task or resting-state conditions. Although our results show that a simple ROI-based correlation analysis may be used to identify significant functional connections in line with the pattern of major descending and ascending projections in humans, future investigations are needed to evaluate the performance of various brain-spinal cord functional connectivity methods.\u003c/p\u003e \u003cp\u003eFurthermore, the improved in-plane spatial resolution at the brain level enabled us to detect significant activation clusters in important brainstem nuclei, including the reticular formation (MRF, and PRF), and red nucleus. Importantly, we detected bilateral activation clusters in PRF and MRF, with stronger activation on the ipsilateral side (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e), consistent with the observed pattern of reticulospinal projections targeting forelimb muscles in primates (Sakai et al., \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Soteropoulos et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). The contribution of various brainstem nuclei in the reticular formation, including the lateral rostral medulla (Arber and Costa, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Ruder et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Ruder and Arber, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), and caudal medulla (Esposito et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) (corresponding to MRF in humans) in forelimb motor control has been recently the topic of intense investigations in rodents. We hope that FASB acquisition/analysis pipeline provides a standard method for examination of task-related activation and functional interaction across cortical, subcortical, brainstem, and spinal cord levels, and enables testing/translation of preclinical findings to humans and clinical populations.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eFunding:\u003c/h2\u003e \u003cp\u003eResearch reported in this study was supported by NIH NINDS under Award Number R01NS131227 (PI: Vahdat), and Florida Department of Health James and Esther King Biomedical Research, Award Number 23K07 (PI: Vahdat). This work was supported in part by an NIH award, S10 OD021726, for High End Instrumentation. JD was supported by the Courtois Foundation, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.canadahelps.org/en/charities/fondation-courtois\u003c/span\u003e\u003cspan address=\"https://www.canadahelps.org/en/charities/fondation-courtois\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. C.L was supported by Fonds de Recherche du Qu\u0026eacute;bec\u0026mdash;Sant\u0026eacute; (FRQ-S) and the Canadian Institutes of Health Research (CIHR).\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eS.V. designed the experiments, optimized the pipeline, and wrote the codes. S.V., F.B., C.L. collected the data, performed data analysis and generated the figures. S.V., O.L., B.D., J.D., F.B., and C.L wrote and revised the manuscript. All authors reviewed the manuscript.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAfyouni, S., Nichols, T.E., 2018. Insight and inference for DVARS. Neuroimage 172, 291\u0026ndash;312. https://doi.org/10.1016/J.NEUROIMAGE.2017.12.098\u003c/li\u003e\n\u003cli\u003eArber, S., Costa, R.M., 2022. Networking brainstem and basal ganglia circuits for movement. Nat. Rev. Neurosci. 23, 342\u0026ndash;360. https://doi.org/10.1038/S41583-022-00581-W\u003c/li\u003e\n\u003cli\u003eArcher, D.B., Kang, N., Misra, G., Marble, S., Patten, C., Coombes, S.A., 2018. Visual feedback alters force control and functional activity in the visuomotor network after stroke. NeuroImage Clin. 17, 505\u0026ndash;517. https://doi.org/10.1016/j.nicl.2017.11.012\u003c/li\u003e\n\u003cli\u003eArcher, D.B., Misra, G., Patten, C., Coombes, S.A., 2016. Microstructural properties of premotor pathways predict visuomotor performance in chronic stroke. Hum. Brain Mapp. 37, 2039\u0026ndash;2054. https://doi.org/10.1002/hbm.23155\u003c/li\u003e\n\u003cli\u003eBaker, S.N., 2011. The primate reticulospinal tract, hand function and functional recovery. J. Physiol. 589, 5603\u0026ndash;5612. https://doi.org/10.1113/JPHYSIOL.2011.215160\u003c/li\u003e\n\u003cli\u003eBarry, R.L., Conrad, B.N., Smith, S.A., Gore, J.C., 2018a. A practical protocol for measurements of spinal cord functional connectivity. Sci. Rep. 8. https://doi.org/10.1038/s41598-018-34841-6\u003c/li\u003e\n\u003cli\u003eBarry, R.L., Vannesjo, S.J., By, S., Gore, J.C., Smith, S.A., 2018b. Spinal cord MRI at 7T. Neuroimage 168, 437\u0026ndash;451. https://doi.org/10.1016/J.NEUROIMAGE.2017.07.003\u003c/li\u003e\n\u003cli\u003eBehzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage 37, 90\u0026ndash;101. https://doi.org/10.1016/J.NEUROIMAGE.2007.04.042\u003c/li\u003e\n\u003cli\u003eBraa\u0026szlig;, H., Feldheim, J., Chu, Y., Tinnermann, A., Finsterbusch, J., B\u0026uuml;chel, C., Schulz, R., Gerloff, C., 2023. Association between activity in the ventral premotor cortex and spinal cord activation during force generation\u0026mdash;A combined cortico‐spinal fMRI study. Hum. Brain Mapp. 44, 6471. https://doi.org/10.1002/HBM.26523\u003c/li\u003e\n\u003cli\u003eBrooks, J.C.W., Beckmann, C.F., Miller, K.L., Wise, R.G., Porro, C.A., Tracey, I., Jenkinson, M., 2008. Physiological noise modelling for spinal functional magnetic resonance imaging studies. Neuroimage 39, 680\u0026ndash;692. https://doi.org/10.1016/j.neuroimage.2007.09.018\u003c/li\u003e\n\u003cli\u003eBurciu, R.G., Chung, J.W., Shukla, P., Ofori, E., Li, H., McFarland, N.R., Okun, M.S., Vaillancourt, D.E., 2016. Functional MRI of disease progression in Parkinson disease and atypical parkinsonian syndromes. Neurology 87, 709\u0026ndash;717. https://doi.org/10.1212/WNL.0000000000002985\u003c/li\u003e\n\u003cli\u003eCauzzo, S., Singh, K., Stauder, M., Garc\u0026iacute;a-Gomar, M.G., Vanello, N., Passino, C., Staab, J., Indovina, I., Bianciardi, M., 2022. Functional connectome of brainstem nuclei involved in autonomic, limbic, pain and sensory processing in living humans from 7 Tesla resting state fMRI. Neuroimage 250. https://doi.org/10.1016/J.NEUROIMAGE.2022.118925\u003c/li\u003e\n\u003cli\u003eCox, R.W., 1996. AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res. 29, 162\u0026ndash;173. https://doi.org/10.1006/cbmr.1996.0014\u003c/li\u003e\n\u003cli\u003eDavidson, A.G., Schieber, M.H., Buford, J.A., 2007. Bilateral spike-triggered average effects in arm and shoulder muscles from the monkey pontomedullary reticular formation. J. Neurosci. 27, 8053\u0026ndash;8058. https://doi.org/10.1523/JNEUROSCI.0040-07.2007\u003c/li\u003e\n\u003cli\u003eDe Leener, B., Fonov, V.S., Collins, D.L., Callot, V., Stikov, N., Cohen-Adad, J., 2018. PAM50: Unbiased multimodal template of the brainstem and spinal cord aligned with the ICBM152 space. Neuroimage 165, 170\u0026ndash;179. https://doi.org/10.1016/J.NEUROIMAGE.2017.10.041\u003c/li\u003e\n\u003cli\u003eDe Leener, B., L\u0026eacute;vy, S., Dupont, S.M., Fonov, V.S., Stikov, N., Louis Collins, D., Callot, V., Cohen-Adad, J., 2017. SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data. Neuroimage 145, 24\u0026ndash;43. https://doi.org/10.1016/j.neuroimage.2016.10.009\u003c/li\u003e\n\u003cli\u003eDoyon, J., Gabitov, E., Vahdat, S., Lungu, O., Boutin, A., 2018. Current issues related to motor sequence learning in humans. Curr. Opin. Behav. Sci. 20, 89\u0026ndash;97. https://doi.org/10.1016/J.COBEHA.2017.11.012\u003c/li\u003e\n\u003cli\u003eEsposito, M.S., Capelli, P., Arber, S., 2014. Brainstem nucleus MdV mediates skilled forelimb motor tasks. Nature 508, 351\u0026ndash;356. https://doi.org/10.1038/nature13023\u003c/li\u003e\n\u003cli\u003eFinsterbusch, J., Eippert, F., Buchel, C., 2012. Single, slice-specific z-shim gradient pulses improve T2*-weighted imaging of the spinal cord. Neuroimage 59, 2307\u0026ndash;2315. https://doi.org/10.1016/j.neuroimage.2011.09.038\u003c/li\u003e\n\u003cli\u003eFinsterbusch, J., Sprenger, C., Buchel, C., 2013. Combined T2*-weighted measurements of the human brain and cervical spinal cord with a dynamic shim update. Neuroimage 79, 153\u0026ndash;161. https://doi.org/10.1016/j.neuroimage.2013.04.021\u003c/li\u003e\n\u003cli\u003eGlover, G.H., Li, T.Q., Ress, D., 2000. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn. Reson. Med. 44, 162\u0026ndash;167. https://doi.org/10.1002/1522-2594(200007)44:1\u0026lt;162::AID-MRM23\u0026gt;3.0.CO;2-E\u003c/li\u003e\n\u003cli\u003eGros, C., De Leener, B., Badji, A., Maranzano, J., Eden, D., Dupont, S.M., Talbott, J., Zhuoquiong, R., Liu, Y., Granberg, T., Ouellette, R., Tachibana, Y., Hori, M., Kamiya, K., Chougar, L., Stawiarz, L., Hillert, J., Bannier, E., Kerbrat, A., Edan, G., Labauge, P., Callot, V., Pelletier, J., Audoin, B., Rasoanandrianina, H., Brisset, J.C., Valsasina, P., Rocca, M.A., Filippi, M., Bakshi, R., Tauhid, S., Prados, F., Yiannakas, M., Kearney, H., Ciccarelli, O., Smith, S., Treaba, C.A., Mainero, C., Lefeuvre, J., Reich, D.S., Nair, G., Auclair, V., McLaren, D.G., Martin, A.R., Fehlings, M.G., Vahdat, S., Khatibi, A., Doyon, J., Shepherd, T., Charlson, E., Narayanan, S., Cohen-Adad, J., 2019. Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks. Neuroimage 184, 901\u0026ndash;915. https://doi.org/10.1016/j.neuroimage.2018.09.081\u003c/li\u003e\n\u003cli\u003eHyvarinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Netw 10, 626\u0026ndash;634. https://doi.org/10.1109/72.761722\u003c/li\u003e\n\u003cli\u003eIslam, H., Law, C.S.W., Weber, K.A., Mackey, S.C., Glover, G.H., 2019. Dynamic per slice shimming for simultaneous brain and spinal cord fMRI. Magn. Reson. Med. 81, 825\u0026ndash;838. https://doi.org/10.1002/MRM.27388\u003c/li\u003e\n\u003cli\u003eKaptan, M., Horn, U., Vannesjo, S.J., Mildner, T., Weiskopf, N., Finsterbusch, J., Brooks, J.C.W., Eippert, F., 2023. Reliability of resting-state functional connectivity in the human spinal cord: Assessing the impact of distinct noise sources. Neuroimage 275. https://doi.org/10.1016/J.NEUROIMAGE.2023.120152\u003c/li\u003e\n\u003cli\u003eKasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. J. Neurosci. Methods 276, 56\u0026ndash;72. https://doi.org/10.1016/J.JNEUMETH.2016.10.019\u003c/li\u003e\n\u003cli\u003eKhatibi, A., Vahdat, S., Lungu, O., Finsterbusch, J., B\u0026uuml;chel, C., Cohen-Adad, J., Marchand-Pauvert, V., Doyon, J., 2022. Brain-spinal cord interaction in long-term motor sequence learning in human: An fMRI study. Neuroimage 253, 119111. https://doi.org/10.1016/J.NEUROIMAGE.2022.119111\u003c/li\u003e\n\u003cli\u003eKinany, N., Khatibi, A., Lungu, O., Finsterbusch, J., B\u0026uuml;chel, C., Marchand-Pauvert, V., Van De Ville, D., Vahdat, S., Doyon, J., 2023. Decoding cerebro-spinal signatures of human behavior: Application to motor sequence learning. Neuroimage 275. https://doi.org/10.1016/J.NEUROIMAGE.2023.120174\u003c/li\u003e\n\u003cli\u003eKinany, N., Pirondini, E., Micera, S., Van De Ville, D., 2022. Spinal Cord fMRI: A New Window into the Central Nervous System. Neuroscientist. https://doi.org/10.1177/10738584221101827\u003c/li\u003e\n\u003cli\u003eKong, Y., Jenkinson, M., Andersson, J., Tracey, I., Brooks, J.C.W., 2012. Assessment of physiological noise modelling methods for functional imaging of the spinal cord. Neuroimage 60, 1538\u0026ndash;1549. https://doi.org/10.1016/j.neuroimage.2011.11.077\u003c/li\u003e\n\u003cli\u003eLandelle, C., Dahlberg, L.S., Lungu, O., Misic, B., De Leener, B., Doyon, J., 2023. Altered Spinal Cord Functional Connectivity Associated with Parkinson\u0026rsquo;s Disease Progression. Mov. Disord. 38, 636\u0026ndash;645. https://doi.org/10.1002/MDS.29354\u003c/li\u003e\n\u003cli\u003eLandelle, C., Lungu, O., Vahdat, S., Kavounoudias, A., Marchand-Pauvert, V., De Leener, B., Doyon, J., 2021. Investigating the human spinal sensorimotor pathways through functional magnetic resonance imaging. Neuroimage 245. https://doi.org/10.1016/J.NEUROIMAGE.2021.118684\u003c/li\u003e\n\u003cli\u003eLemon, R.N., 2008. Descending pathways in motor control. Annu Rev Neurosci 31, 195\u0026ndash;218. https://doi.org/10.1146/annurev.neuro.31.060407.125547\u003c/li\u003e\n\u003cli\u003eMarrelec, G., Krainik, A., Duffau, H., Pelegrini-Issac, M., Lehericy, S., Doyon, J., Benali, H., 2006. Partial correlation for functional brain interactivity investigation in functional MRI. Neuroimage 32, 228\u0026ndash;237. https://doi.org/10.1016/j.neuroimage.2005.12.057\u003c/li\u003e\n\u003cli\u003ePierrot-Deseilligny, E., Burke, D., 2012. The circuitry of the human spinal cord: Spinal and corticospinal mechanisms of movement, The Circuitry of the Human Spinal Cord: Spinal and Corticospinal Mechanisms of Movement. https://doi.org/10.1017/CBO9781139026727\u003c/li\u003e\n\u003cli\u003eRaudino, F., Leva, S., 2012. Involvement of the spinal cord in Parkinson\u0026rsquo;s disease. Int. J. Neurosci. 122, 1\u0026ndash;8. https://doi.org/10.3109/00207454.2011.613551\u003c/li\u003e\n\u003cli\u003eRissman, J., Gazzaley, A., D\u0026rsquo;Esposito, M., 2004. Measuring functional connectivity during distinct stages of a cognitive task. Neuroimage 23, 752\u0026ndash;763. https://doi.org/10.1016/J.NEUROIMAGE.2004.06.035\u003c/li\u003e\n\u003cli\u003eRuder, L., Arber, S., 2019. Brainstem Circuits Controlling Action Diversification. Annu. Rev. Neurosci. 42, 485\u0026ndash;504. https://doi.org/10.1146/ANNUREV-NEURO-070918-050201\u003c/li\u003e\n\u003cli\u003eRuder, L., Schina, R., Kanodia, H., Valencia-Garcia, S., Pivetta, C., Arber, S., 2021. A functional map for diverse forelimb actions within brainstem circuitry. Nature 590, 445\u0026ndash;450. https://doi.org/10.1038/S41586-020-03080-Z\u003c/li\u003e\n\u003cli\u003eSakai, S.T., Davidson, A.G., Buford, J.A., 2009. Reticulospinal neurons in the pontomedullary reticular formation of the monkey (Macaca fascicularis). Neuroscience 163, 1158\u0026ndash;1170. https://doi.org/10.1016/J.NEUROSCIENCE.2009.07.036\u003c/li\u003e\n\u003cli\u003eSingh, K., Garc\u0026iacute;a-Gomar, M.G., Bianciardi, M., 2021. Probabilistic Atlas of the Mesencephalic Reticular Formation, Isthmic Reticular Formation, Microcellular Tegmental Nucleus, Ventral Tegmental Area Nucleus Complex, and Caudal-Rostral Linear Raphe Nucleus Complex in Living Humans from 7 Tesla Magnetic Resonance Imaging. Brain Connect. 11, 613\u0026ndash;623. https://doi.org/10.1089/BRAIN.2020.0975\u003c/li\u003e\n\u003cli\u003eSmith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E., Johansen-Berg, H., Bannister, P.R., De Luca, M., Drobnjak, I., Flitney, D.E., Niazy, R.K., Saunders, J., Vickers, J., Zhang, Y., De Stefano, N., Brady, J.M., Matthews, P.M., 2004. Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23 Suppl 1, S208-19. https://doi.org/S1053-8119(04)00393-3 [pii]10.1016/j.neuroimage.2004.07.051\u003c/li\u003e\n\u003cli\u003eSoteropoulos, D.S., Williams, E.R., Baker, S.N., 2012. Cells in the monkey ponto-medullary reticular formation modulate their activity with slow finger movements. J. Physiol. 590, 4011\u0026ndash;4027. https://doi.org/10.1113/jphysiol.2011.225169\u003c/li\u003e\n\u003cli\u003eTinnermann, A., B\u0026uuml;chel, C., Cohen-Adad, J., 2021. Cortico-spinal imaging to study pain. Neuroimage. https://doi.org/10.1016/j.neuroimage.2020.117439\u003c/li\u003e\n\u003cli\u003eVahdat, S., Darainy, M., Milner, T.E., Ostry, D.J., 2011. Functionally specific changes in resting-state sensorimotor networks after motor learning. J Neurosci 31, 16907\u0026ndash;16915. https://doi.org/10.1523/JNEUROSCI.2737-11.2011\u003c/li\u003e\n\u003cli\u003eVahdat, S., Darainy, M., Ostry, D.J., 2014. Structure of Plasticity in Human Sensory and Motor Networks Due to Perceptual Learning. J. Neurosci. 34, 2451\u0026ndash;2463. https://doi.org/10.1523/JNEUROSCI.4291-13.2014\u003c/li\u003e\n\u003cli\u003eVahdat, S., Fogel, S., Benali, H., Doyon, J., 2017. Network-wide reorganization of procedural memory during NREM sleep revealed by fMRI. Elife 6, e24987. https://doi.org/10.7554/eLife.24987\u003c/li\u003e\n\u003cli\u003eVahdat, S., Khatibi, A., Lungu, O., Finsterbusch, J., B\u0026uuml;chel, C., Cohen-Adad, J., Marchand-Pauvert, V., Doyon, J., 2020. Resting-state brain and spinal cord networks in humans are functionally integrated. PLoS Biol. 18(7): e30. https://doi.org/10.1371/journal.pbio.3000789\u003c/li\u003e\n\u003cli\u003eVahdat, S., Lungu, O., Cohen-Adad, J., Marchand-Pauvert, V., Benali, H., Doyon, J., 2015. Simultaneous brain\u0026ndash;cervical cord fMRI reveals intrinsic spinal cord plasticity during motor sequence learning. PLoS Biol. 13, 1\u0026ndash;25. https://doi.org/10.1371/journal.pbio.1002186\u003c/li\u003e\n\u003cli\u003eVahdat, S., Maneshi, M., Grova, C., Gotman, J., Milner, T.E., 2012. Shared and specific independent components analysis for between-group comparison. Neural Comput. 24, 3052\u0026ndash;3090. https://doi.org/10.1162/NECO_a_00355\u003c/li\u003e\n\u003cli\u003eVaillancourt, D.E., Thulborn, K.R., Corcos, D.M., 2003. Neural Basis for the Processes That Underlie Visually Guided and Internally Guided Force Control in Humans. J. Neurophysiol. 90, 3330\u0026ndash;3340. https://doi.org/10.1152/jn.00394.2003\u003c/li\u003e\n\u003cli\u003eWahl, A.S., Schwab, M.E., 2014. Finding an optimal rehabilitation paradigm after stroke: Enhancing fiber growth and training of the brain at the right moment. Front. Hum. Neurosci. 8, 381. https://doi.org/10.3389/FNHUM.2014.00381/BIBTEX\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-3889284/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3889284/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eSimultaneous functional magnetic resonance imaging (fMRI) of the spinal cord and brain represents a powerful method for examining both ascending sensory and descending motor pathways in humans \u003cem\u003ein vivo\u003c/em\u003e. However, its image acquisition protocols, and processing pipeline are less well established. This limitation is mainly due to technical difficulties related to spinal cord fMRI, and problems with the logistics stemming from a large field of view covering both brain and cervical cord. Here, we propose an acquisition protocol optimized for both anatomical and functional images, as well as an optimized integrated image processing pipeline, which consists of a novel approach for automatic modeling and mitigating the negative impact of spinal voxels with low temporal signal to noise ratio (tSNR). We validate our integrated pipeline, named FASB, using simultaneous fMRI data acquired during the performance of a motor task, as well as during resting-state conditions. We demonstrate that FASB outperforms the current spinal fMRI processing methods in three domains, including motion correction, registration to the spinal cord template, and improved detection power of the group-level analysis by removing the effects of participant-specific low tSNR voxels, typically observed at the disk level. Using FASB, we identify significant task-based activations in the expected sensorimotor network associated with a unilateral handgrip force production task across the entire central nervous system, including the contralateral sensorimotor cortex, thalamus, striatum, cerebellum, brainstem, as well as ipsilateral ventral horn at C5-C8 cervical levels. Additionally, our results show significant task-based functional connectivity between the key sensory and motor brain areas and the dorsal and ventral horns of the cervical cord. Overall, our proposed acquisition protocol and processing pipeline provide a robust method for characterizing the activation and functional connectivity of distinct cortical, subcortical, brainstem and spinal cord regions in humans.\u003c/p\u003e","manuscriptTitle":"FASB: an integrated processing pipeline for Functional Analysis of simultaneous Spinal cord-Brain fMRI","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-01-30 18:24:08","doi":"10.21203/rs.3.rs-3889284/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"a7665f6e-3532-4174-acbb-09fcfde976d4","owner":[],"postedDate":"January 30th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":28432904,"name":"Biological sciences/Biological techniques/Imaging/Functional magnetic resonance imaging"},{"id":28432905,"name":"Biological sciences/Neuroscience/Neural circuit"},{"id":28432906,"name":"Biological sciences/Neuroscience/Sensorimotor processing"}],"tags":[],"updatedAt":"2024-02-20T09:45:13+00:00","versionOfRecord":[],"versionCreatedAt":"2024-01-30 18:24:08","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-3889284","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3889284","identity":"rs-3889284","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}