Measurement of stretch-evoked brainstem function using fMRI

Knowledge on the organization of motor function in the reticulospinal tract (RST) is limited by the lack of methods for measuring RST function in humans. Behavioral studies suggest the involvement of the RST in long latency responses (LLRs). LLRs, elicited by precisely controlled perturbations, can therefore act as a viable paradigm to measure motor-related RST activity using functional Magnetic Resonance Imaging (fMRI). Here we present StretchfMRI, a novel technique developed to study RST function associated with LLRs. StretchfMRI combines robotic perturbations with electromyography and fMRI to simultaneously quantify muscular and neural activity during stretch-evoked LLRs without loss of reliability. Using StretchfMRI, we established the muscle-specific organization of LLR activity in the brainstem. The observed organization is partially consistent with animal models, with activity primarily in the ipsilateral medulla for flexors and in the contralateral pons for extensors, but also includes other areas, such as the midbrain and bilateral pontomedullary contributions.


Supplementary analysis of fMRI data
The GLM selected to partition the variance of the measured BOLD signal (Eq. 4 of the main manuscript) does not include separate regressors for background activity and LLR amplitude. Because participants are asked to contract their muscles prior to each perturbation, there is an intrinsic statistical association between the muscle state before the perturbation (i.e., background muscle activity) and the perturbation itself (i.e., reflex response). As such, given the poor temporal resolution of fMRI, it is possible that some of the variance in the BOLD signal explained by LLR-related regressors might actually arise from the neural activity associated to background contraction. To rule out this concern, we conducted two analyses, one based on simulated data, and one based on the implementation of a different GLM on the fMRI data collected in Experiment 2.

Model-based analysis
We initially sought to quantify whether the implemented analysis method is sufficiently sensitive and specific in identifying neural activity related to LLR in presence of correlated background activity. To this goal we have conducted a simulated signal analysis following the methods presented in our previous work 1 .

Methods
We simulated nine different representative brain voxels V 1 to V 9 that selectively respond only to linear combinations of conditions of flexor and extensor muscle activity, defined in terms of LLR events (H F and H E , respectively for flexor and extensor) and background contraction (B F and B E , respectively for flexor and extensor). Specifically, H F and H E are obtained by convolution of rectangular functions (duration: 50 ms, onset: 50 ms after perturbation start, amplitude:H FCR | st ANC,P4 andH ECU | st ANC,P4 , respectively) with the standard HRF (Fig. S8). Similarly, the terms B F and B E are given by the convolution of rectangular functions (unitary magnitude in the time interval when the measured torque has a magnitude greater than 25% of the desired torque) with the standard HRF (Fig. S8). To simulate all four conditions, we used the EMG and torque data measured in a representative participant (participant P4). A schematic representation of the simulated voxels is reported in Tabs. S1 and S2. Table S1. Definition of voxels used in the simulation analysis

Effects
Flexor Extensor Both For each of the simulated voxels, we then added zero-mean Gaussian noise and down-sampled the signal with a period T R = 1.225 s, to simulate realistic BOLD measurements specific to either cortical voxels (Signal-To-Noise ratio equal to 1.2) or brainstem voxels (Signal-to-Noise ratio equal to 0.5). Finally, we quantified the sensitivity that the two proposed GLMs (GLM 1 in Eq. 4 and GLM 2 in Eq. S1) have in identifying the four different regressors, by computing the true positive (TPR) and false positive (FPR) rates afforded when either of the two GLM is used to fit the signal simulated in each of the nine different voxels.

Results
The results show that the model GLM 1 affords the highest TPR for LLR-specific voxels in both cortical and brainstem voxels, and for both flexor and extensor regressors (V 1 and V 2 , respectively; Fig. S3A, top row and Table S2). However, it also returns high FPR for LLR-related activity in voxels whose activity is associated to the application fo background torque (Fig. S3A,  middle and Table S2), though the FPR is smaller in the brainstem voxels than in the cortical ones.
The inclusion of the background-specific regressors drastically reduces the FPR of LLR-related activity in backgroundspecific voxels for both flexor and extensor (V 4 and V 5 , respectively; Fig. S3B, middle row and Table S2). However, while TPR remains high in LLR-specific voxels (Fig. S3A, top row and Table S2), it decreases to almost zero in voxels that responds to both LLR and background activity (V 7 and V 8 ; Fig. S3A, bottom row and Table S2).

Alternative GLM analysis
We conducted an alternative statistical analysis for the BOLD fMRI data that aims to determine the LLR-related neural activity including a correction for the possible confounding effect introduced by the background-related activity. Figure S2. Group-average of the flexion/extension angle (A) and torque (B) measured in response to RaH perturbations applied at different velocities. In all graph the thick line represents the mean and the shaded area the 95% confidence interval. To facilitate interpretation on the graphs a gray shaded area has been included in each graph in the time interval where a LLR is expected. Figure S3. Identification rate for GLM 1 (A) and GLM 2 (B). Bold labels indicate true positive association between the indicated regressor and the BOLD signal measured in the corresponding voxel. Table S2. Percent state detected in the different voxels for the two GLMs. Bolded labels represents true positive identification rates, while un-bolded labels represents false positive identification rate .

Voxel
Effects Bolded number represents true positive identification rates, while un-bolded numbers represents false positive identification rates

Methods
Methods closely followed those described in Sec. 2.5 for most aspects. Separately for each session, we performed two first-level analyses to determine activation in the whole brain and in the brainstem. While we used the same general linear model (GLM) for the whole-brain and brainstem analyses, we used regionally-specific hemodynamic response functions (HRF) for different brain regions. A brainstem-specific HRF 2 was used to assess activation in the brainstem, while a standard HRF was used to assess activation in all other regions. In this alternative GLM (GLM 2 ) analysis, the neural response was modeled as: where H FCR and H ECU are obtained by convolution of rectangular functions (duration: 50 ms, onset: 50 ms after perturbation start, amplitude:H FCR | st ANC andH ECU | st ANC , respectively) with the appropriate HRF (Fig. S8). The terms B FCR and B ECU are now added to the model to account for metabolic activity associated with muscle-specific pre-perturbation (i.e., background) contraction state. They are similarly obtained convolving rectangular functions (unitary magnitude in the time interval when the measured torque is in magnitude greater than 25% of the desired torque) with the appropriate HRF (Fig. S8). As in all analyses, a set of 6 nuisance regressors R was included to account for variance in the measured signal associated with the 3D head movements (translation and rotations).

Results
For the brainstem-specific analysis, muscle-specific maps included several clusters of activation associated only to the LLR-specific regressors (Fig. S4, Table S3). In fact, no supra-threshold voxels were identified to correlate with the backgroundspecific regressors (Fig. S5). LLR-related activity specific to FCR included a total of 18 voxels, including two clusters spanning bilaterally the midbrain, two clusters spanning bilaterally the pons, and two clusters in the right superior medulla (Fig. S4 top, and Table S3). LLR-related activity specific to ECU included a total of 16 voxels, including one voxel in the right midbrain, four clusters spanning bilaterally the pons, and one cluster in the right medulla (Fig. S4 top, and Table S3).
In contrast to what we observed in the brainstem-specific analysis, in the whole brain analysis muscle-specific maps included several clusters of activation associated to both LLR-specific ( Fig. S4 bottom, Table S4, and Table S5) and background-specific regressors ( Fig. S4 bottom, Table S6, and Table S7). Specifically, cortical activation associated to the two background-specific regressors (B FCR and B ECU ) primarily included overlapped clusters of activation in the contralateral sensorimotor cortex and 4/30 inferior parietal lobule, in the ipsilateral insular and opercular cortex, and in the bilateral visual cortex. Cerebellar activity included overlapped clusters in the ipsilateral region V and contralateral region VI, while in the subcortical regions activity was primarily detected in the contralateral putamen ( Fig. S4 bottom, Table S6, and Table S7). LLR-related cortical activity specific to FCR primarily included contralateral activation in the sensorimotor cortex, inferior parietal lobule and bilateral activation in the visual cortex. Cerebellar activity could be instead observed bilaterally in the regions I-IV and VIII, and contralaterally in the region VI and in the Crus II, while subcortical activity could be observed in the contralateral putamen and thalamus and in the ipsilateral caudate nucleus (Fig. S4 bottom, Table S4). As for the LLR-related activity specific to ECU, no subcortical activity could be observed, with cortical activity observable only in the ipsilateral visual cortex. Cerebellar activity was instead more diffuse and included bilateral activation in the region VIII and Crus I, ipsilateral activity in the regions V and X, and contralateral activity in the region VI and in the Crus II ( Fig. S4 bottom, Table S4). Figure S4. (Top) LLR-specific activation maps in the brainstem for FCR and ECU. The maps refer to contrasts β FCR,L > 0 and β ECU,L > 0 obtained for the brainstem-specific analysis. The threshold t-statistic, obtained after FWE correction, was equal to 5.62 for both regressors. For reference, the t-statistic that refer to a Bonferroni correction (t Bon = 5.8) is marked on each colorbar with a black line. (Bottom) LLR-specific activation maps in the whole-brain for FCR and ECU. The threshold t-statistic, obtained after FWE correction, was set to 7.28 and 6.99, respectively for FCR and ECU, while t Max was equal to 12.81 and 11.95, respectively. Colorbars are saturated at t = 10 for better visualization of t-statistic gradients. All statistical parametric maps are overlaid on axial slices of the standard Montreal Neurological Institute 152 template, with reported z coordinate in mm 5/30 Figure S5. (Top) Activation maps in the brainstem for FCR and ECU background activity. The maps refer to contrasts β FCR,B > 0 and β ECU,B > 0 obtained for the brainstem-specific analysis. The threshold t-statistic, obtained after FWE correction, was equal to 5.62 for both regressors. For reference, the t-statistic that refer to a Bonferroni correction (t Bon = 5.8) is marked on each colorbar with a black line. (Bottom) Activation maps in the whole brain for FCR and ECU background activity. The threshold t-statistic, obtained after FWE correction, was set to 7.27 and 7.17, respectively for FCR and ECU, while t Max was equal to 11.44 and 16.66, respectively. Colorbars are saturated at t = 10 for better visualization of t-statistic gradients. All statistical parametric maps are overlaid on axial slices of the standard Montreal Neurological Institute 152 template, with reported z coordinate in mm

Discussion
The results of the simulation analysis show that voxels that are positively associated with the regressor H F in GLM 1 might respond to LLR activity, background activity, or both-i.e. V 1 , V 4 , or V 7 , respectively (Fig. S3A, first column). On the contrary, when using GLM 2 , voxels that are positively associated with the regressor H F can only respond to LLR-related activity-i.e. V 1 -since a 0% positive identification rate characterize voxels V 4 and V 7 . For all three voxels, however, the background-specific regressor B F dominates, reducing the specificity of the regressor B F in voxels that only respond to LLR activity (V 1 ) and the sensitivity of the regressor H F in voxels that instead respond to both LLR and background activity (V 7 ), for both cortical and brainstem voxels.
The same observation can be made considering voxels that respond to activity in the extensors. Therefore, it is possible to conclude that when using the model GLM 2 , background related activity will always be identified if present, and dominate over the LLR-related regressor. However, the brainstem activation maps obtained using the GLM 2 model show the complete absence of voxels whose signal is significantly associated with the background specific regressors (Fig. S5). This evidence, together with the observation that the clusters of activations included in the LLR-specific regressors (Fig. S4, Table S3) are closely aligned with the clusters of activations identified by GLM 1 (Fig. 3, Table 1), suggest that the identified brainstem activity is truly associated with LLR-related activity and not with background activation. The insight from the two analyses presented above advances that it is unlikely that some of the variance in the BOLD signal explained by LLR-related regressors might arise from the neural activity associated to background contraction.
In contrast to the similarity observed between the results obtained for the two GLMs in the brainstem-specific analysis, at the whole brain level the two GLMs produce fairly different results. Since both LLRs and isometric contractions are intrinsic motor actions, they are expected to share most of the same cortical substrates. As such, as predicted by the model, it is possible to observe a decrease of the sensitivity that the model has in identifying LLR-related activity, mostly evident in ECU-specific maps (Fig. S4 bottom, Table S4).
In summary, due to the co-linearity between the LLR-specific and background-specific regressors-i.e. after the application of a background torque there is always a perturbation that elicits a reflex response in the pre-conditioned muscles-our experimental protocol does not afford enough statistical power to reliably decouple the LLR-related and background-related activity at the cortical level. However, as demonstrated by the combination of the results of the simulation analysis and of the alternative GLM analysis, it is very likely that the measured activations are specific to LLR-related responses, and minimally affected by statistically associated activity arising from background contraction.  Figure S6. Stretch-evoked EMG signal measured using off-the-shelf electrodes, amplifier, and signal processing pipeline in response to a 125 deg/s RaH perturbation for the Flexor Carpi Radialis (FCR) for a representative subject.  Plots are obtained with a combination of the novel data collection hardware with a standard signal processing pipeline (left), a time-domain subtraction of measured and reference signals (center), and the novel adaptive noise cancellation pipeline (right). In all graph the thick line represents the mean and the shaded area the 95% confidence interval. To facilitate interpretation on the graphs a gray shaded area has been included in each graph in the time interval where a LLR is expected.       Figure S13. Average BOLD signal measured in response to muscle stretch in the selected voxels, for the selected participants. Specifically, A include the data that refer to the subjects that showed the greatest difference between the statistical scores specific to the FCR and ECU regressors, while B include the data that refer to those subject that showed similar t score for the two regressors (Fig. 6). In all graphs, the solid line shows the average residual BOLD response after variance associated to nuisance regressors is removed, while the shaded area represents the 95% confidence interval of the mean response. Dashed lines show the average modeled response respectively for the LLR-specific regressors (green or black line), and for the background-specific regressor (light blue line). Modeled response for the LLR-specific regressor was scaled based on the appropriate β coefficient estimated by the GLM model, while modeled background response, which is not included in the primary GLM, was arbitrarily scaled to match the same amplitude of the LLR-specific regressor. Measured and LLR-specific signals are then color-coded to distinguish the response that is significant at the group level for that specific muscle (green), from the one that is not significant (black).