Effects of Perturbation Velocity, Direction, Background Muscle Activation, and Task Instruction on Long-Latency Responses Measured From Forearm Muscles

The central nervous system uses feedback processes that occur at multiple time scales to control interactions with the environment. The long-latency response (LLR) is the fastest process that directly involves cortical areas, with a motoneuron response measurable 50 ms following an imposed limb displacement. Several behavioral factors concerning perturbation mechanics and the active role of muscles prior or during the perturbation can modulate the long-latency response amplitude (LLRa) in the upper limbs, but the interactions among many of these factors had not been systematically studied before. We conducted a behavioral study on thirteen healthy individuals to determine the effect and interaction of four behavioral factors – background muscle torque, perturbation direction, perturbation velocity, and task instruction – on the LLRa evoked from the flexor carpi radialis (FCR) and extensor carpi ulnaris (ECU) muscles after velocity-controlled wrist displacements. The effects of the four factors were quantified using both a 0D statistical analysis on the average perturbation-evoked EMG signal in the period corresponding to an LLR, and using a timeseries analysis of EMG signals. All factors significantly modulated LLRa, and their combination nonlinearly contributed to modulating the LLRa. Specifically, all the three-way interaction terms that could be computed without including the interaction between instruction and velocity significantly modulated the LLR. Analysis of the three-way interaction terms of the 0D model indicated that for the ECU muscle, the LLRa evoked when subjects are asked to maintain their muscle activation in response to the perturbations was greater than the one observed when subjects yielded to the perturbations (p < 0.001), but this effect was not measured for muscles undergoing shortening or in absence of background muscle activation. Moreover, higher perturbation velocity increased the LLRa evoked from the stretched muscle in presence of a background torque (p < 0.001), but no effects of velocity were measured in absence of background torque. Also, our analysis identified significant modulations of LLRa in muscles shortened by the perturbation, including an interaction between torque and velocity, and an effect of both torque and velocity. The time-series analysis indicated the significance of additional transient effects in the LLR region for muscles undergoing shortening.

The central nervous system uses feedback processes that occur at multiple time scales to control interactions with the environment. The long-latency response (LLR) is the fastest process that directly involves cortical areas, with a motoneuron response measurable 50 ms following an imposed limb displacement. Several behavioral factors concerning perturbation mechanics and the active role of muscles prior or during the perturbation can modulate the long-latency response amplitude (LLRa) in the upper limbs, but the interactions among many of these factors had not been systematically studied before. We conducted a behavioral study on thirteen healthy individuals to determine the effect and interaction of four behavioral factors -background muscle torque, perturbation direction, perturbation velocity, and task instruction -on the LLRa evoked from the flexor carpi radialis (FCR) and extensor carpi ulnaris (ECU) muscles after velocity-controlled wrist displacements. The effects of the four factors were quantified using both a 0D statistical analysis on the average perturbation-evoked EMG signal in the period corresponding to an LLR, and using a timeseries analysis of EMG signals. All factors significantly modulated LLRa, and their combination nonlinearly contributed to modulating the LLRa. Specifically, all the three-way interaction terms that could be computed without including the interaction between instruction and velocity significantly modulated the LLR. Analysis of the three-way interaction terms of the 0D model indicated that for the ECU muscle, the LLRa evoked when subjects are asked to maintain their muscle activation in response to the perturbations was greater than the one observed when subjects yielded to the perturbations (p < 0.001), but this effect was not measured for muscles undergoing shortening or in absence of background muscle activation. Moreover, higher perturbation velocity increased the LLRa evoked from the stretched muscle in presence of a background torque (p < 0.001), but no effects of velocity were measured in absence of background torque. Also, our analysis identified

INTRODUCTION
Countering unexpected and unpredictable loads is a ubiquitous occurrence of everyday life. Humans can precisely perform movements and interact with the environment even in the presence of these external perturbations. These mechanical perturbations require the nervous system to induce a compensatory action in order to ensure the task success. An important component of the compensatory actions produced by the central nervous system is the long-latency response (LLR). In upper limb muscles, the LLR is evident as the burst of muscle activity occurring 50-100 ms following a limb displacement. Accordingly, this event occurs between the fastest nervous system response, i.e., the short-latency reflex (SLR) occurring within 20-50 ms, and the delayed voluntary reaction which begins 100 ms after the imposed perturbation (Hammond, 1956;Lee and Tatton, 1975;Kurtzer, 2015).
After a seminal study by Hammond in 1956(Hammond, 1956, several investigators have utilized a limb perturbation paradigm to investigate the physiological mechanisms subserving the muscle stretch responses to the externally applied loads (Allum, 1975;Crago et al., 1976;Evarts and Granit, 1976;Thomas et al., 1977;Kurtzer, 2015;Zonnino et al., 2019). In these paradigms, the muscle responses including the LLRs are recorded through surface electromyography (EMG) activity evoked in the muscle stretched by an imposed angular joint displacement induced by a mechanical perturbation of known and controllable characteristics (Rothwell et al., 1980;Tarkka and Larsen, 1987;Cody and Plant, 1989;Matthews, 1989Matthews, , 1993Noth et al., 1991;Kurtzer, 2015).
There is a body of evidence supporting the practical importance of muscle stretch responses -specifically the LLR component -in the neurological research. Previous studies showed that LLR can be considered as the primary outcome measure in various rehabilitation and robot-aided training protocols for several neurological diseases including stroke, Parkinson's disease (PD), spinal cord injury, and cerebellar ataxia (Sinkjaer and Hayashi, 1989;Hayashi et al., 2001;Trumbower et al., 2013;Mirbagheri et al., 2015;Banks et al., 2019;Deneri et al., 2020). One recent study in 2019 (Banks et al., 2019), distinguished LLR as a promising physiological marker of walking dysfunction in chronic stroke. Trumbower et al., also demonstrated that there is a bilateral impaired regulation of the LLR during tasks which require increased stability in both the paretic and nonparetic upper limbs of stroke survivors (Trumbower et al., 2013). Moreover, previous studies on Parkinson's disease reported that LLR might be a useful objective physiological measure of muscle stiffness and rigidity in PD patients (Sinkjaer and Hayashi, 1989;Hayashi et al., 2001).
Several behavioral factors are known to affect the amplitude of LLRs, including the neuromechanical state of the muscle prior perturbation (i.e., muscle length and activation) (Bedingham and Tatton, 1984;Calancie and Bawa, 1985), the direction of perturbation (i.e., whether the perturbation stretches or shortens the muscle) (Miscio et al., 2001;Lewis et al., 2004Lewis et al., , 2010, the kinematic features of the applied perturbation (i.e., perturbation velocity, duration, amplitude, velocity profile) (Lee and Tatton, 1982;Lewis et al., 2005;Schuurmans et al., 2009), and the instructions provided to participants as to how to respond to the applied perturbations (Miscio et al., 2001;Lewis et al., 2006;Kurtzer et al., 2014).
Although investigators mostly focused on studying LLR features in stretched muscles, there is evidence of EMG activity evoked in the muscle shortened by the applied perturbation (Miscio et al., 2001;Lewis et al., 2004). Specifically, an increase in the EMG activity of the extensor carpi radialis (ECR) muscle which was shortened due to the applied wrist extension perturbation was documented in a study conducted in 2004 (Lewis et al., 2004). Although the shortened muscle response was smaller in amplitude than the one evoked in the stretched (flexor) muscle, both had two separate components of SLR and LLR with a similar onset timing. Another study also observed a similarly timed, low-amplitude EMG response with an onset of about 50 ms, evoked in the ECR muscle in response to a rapid wrist extension (Miscio et al., 2001). The authors suggest that part of the measured effects may be due to crosstalk (a volume-conducted response from the stretched muscle). However, because significant perturbation-evoked responses were measured in muscles undergoing shortening even using intramuscular EMG (Lewis et al., 2010), it may be reasonable that LLR are evoked in muscles subject to both a shortening and a stretching perturbation.
Studies also examined the effects of background muscle activation prior to the imposed perturbation on the LLR amplitude. In general, an increase in background activation results in an increase in the magnitude of muscle activity in both the proximal and distal muscles of the upper limbs (Bedingham and Tatton, 1984;Miscio et al., 2001;Pruszynski et al., 2009). Affecting the background motoneuron pool excitability, preexisting background muscle activation is thought to reflect an automatic adjustment mechanism, known as the automatic gain component of the LLR (Bedingham and Tatton, 1984;Matthews, 1986;Miscio et al., 2001;Pruszynski et al., 2009).
Several studies quantified the effects of the kinematic features of the applied perturbation on the LLR amplitude Tatton, 1975, 1982;Lewis et al., 2005Lewis et al., , 2006Schuurmans et al., 2009). It is generally accepted that the LLR amplitude increases as a function of the velocity of the applied perturbations (Tatton and Bawa, 1979;Bedingham and Tatton, 1984). In the common ramp-and-hold perturbation paradigms, which are conducted at constant velocity, perturbation duration may also play a factor when the duration of the perturbation is within the range of neuromuscular delays expected for the LLR. However, the details of the interaction between perturbation velocity and duration in modulating LLR amplitude are not yet completely understood. A study by Lewis et al. showed that LLR amplitudes of the biceps brachii undergoing stretch are modulated by velocity in all conditions, but the slope of the relationship is also modulated by duration (Lewis et al., 2005). Yet, the range of velocities and durations that modulate the response in such a way is likely limited. In fact, we know that for FCR, a very high velocity and short duration (<40 ms) perturbation is not sufficient to evoke a long-latency response, whereas a perturbation of low velocity and long duration (>60 ms) generates well-developed LLRs (Lee and Tatton, 1982).
Task instruction also plays a key role in modulation of the LLR response. Accumulating evidence shows that the temporal overlap of two different responses including a task-dependent response and an automatic response results in the task-dependent change in LLR amplitude (Rothwell et al., 1980;Lewis et al., 2006;Pruszynski et al., 2011). The task-dependent response is larger when participants attempt to counter a perturbation than yield to the perturbation (Calancie and Bawa, 1985;Miscio et al., 2001;Lewis et al., 2006;Kurtzer et al., 2014). Participants can be instructed to respond to the perturbation in different ways: they can be asked to relax immediately following the perturbation (Miscio et al., 2001) -a condition referred to as "yield"; to maintain the background torque and avoid a voluntary response to the perturbation (Calancie and Bawa, 1985) -a condition referred to as "Do Not Intervene"; or to explicitly compensate by activating their muscles in the opposite direction of the perturbation (Lewis et al., 2006), or to compensate cued by a visual feedback of the hand position (Kurtzer et al., 2014) -conditions referred to as "Resist." In conditions where the subject was instructed to counter the stretch, LLR amplitude typically increased compared to a control condition, demonstrating that the LLR can be modulated to functionally adapt to the task at the upper limb (Crago et al., 1976;Colebatch et al., 1979;Rothwell et al., 1980). However, the 'yield' or the 'DNI' instructions were never directly compared, with the exception of one study (Calancie and Bawa, 1985), which only recruited two individuals. Because the state of the muscle in the "yield" and the "DNI" conditions is fundamentally different, it would be useful to quantify the differential effects of these two conditions on the evoked LLR.
Gathered together, several factors concerning the mechanics of the applied perturbations and the active role of muscles prior or during the perturbation can modulate the amplitude of longlatency responses in the upper limbs. Hence, it is of a paramount importance to study how the interaction of these factors would affect the muscle stretch responses in a single study. However, the majority of previous studies has systematically studied only one or two of the factors modulating LLRa, with the interaction between perturbation velocity and task instruction studied in Lewis et al. (2006), the interaction between perturbation velocity and background torque studied in Bedingham and Tatton (1984), and the interaction between task instruction and background torque studied in Calancie and Bawa (1985). One previous study (Miscio et al., 2001) has studied the effects of three of the factors highlighted above (i.e., task instructions, perturbation direction, and torque), though not with a full factorial design capable of quantifying the interactions among all factors. As such, to the best of our knowledge, no previous study conducted a full factorial design capable of quantifying the effects of and the interactions among all combinations of four factors known to modulate LLR amplitude, i.e., task instructions, perturbation duration, background torque, and perturbation velocity.
The goal of this study is to determine the effects of and the interactions among several experimental factors modulating the LLR amplitude during ramp-and-hold perturbation. Specifically, the goal of this study is to establish the effect and interaction of background muscle torque, perturbation direction, perturbation velocity, and task instruction on the LLR amplitude evoked from the flexor carpi radialis (FCR) and extensor carpi ulnaris (ECU) muscles following the application of controlled angular displacements of the wrist in both the extension and flexion directions.

MATERIALS AND METHODS
Thirteen healthy individuals were recruited to participate in this study (protocol approved by the University of Delaware Institutional Review Board, protocol no. 1097082-6). Subjectsage (mean ± s.d.: 24 ± 3 years) were naïve to the purpose of the study and free from known neurological or orthopedic disorders affecting arm function. Subjects were exposed to an experiment that aimed to quantify the amplitude of long-latency responses via the recording of EMG activity from a wrist flexor and extensor pair, the flexor carpi radialis (FCR) and extensor carpi ulnaris (ECU). These responses were evoked by flexion or extension perturbations applied by a robot to a subject's wrist in various conditions.

Materials
The equipment used for this experiment is shown in Figure 1, and includes several components, described below in detail.

Perturbation Robot
A custom-developed robot, the MR-StretchWrist, was used to apply perturbations to subject's wrists. The MR-StretchWrist is a 1-degree of freedom robot that can provide wrist flexion and extension between −45 to 45 degrees (Zonnino et al., 2019). The robot employs an ultrasonic piezoelectric motor (EN6060, Shinsei Motor Inc., Japan) that can provide 500 mNm of torque and can move at velocities of up to 900 degrees/second.
To provide torque for sufficient muscle stretch within the desired time of 50 ms, a capstan transmission with a 3:1 gear ratio was included in the design. The capstan drive contains two pulleys with different diameters connected via a microfiber braided line (SpiderWire Stealth SPW-0039, 0.4 mm diameter braided fishing line). The cable is wrapped around each pulley multiple times to ensure zero slippage. The capstan transmission is an ideal candidate for this application because it has no backlash, low friction, and high bandwidth.
Measurement of the wrist flexion/extension angle was obtained using an incremental encoder (resolution: 0.09 deg) placed on the motor shaft, with a resulting resolution in measuring the wrist flexion/extension angle of 0.03 deg.

EMG Amplifier and Electrodes
Electromyography data was recorded with an OTBioelettronica EMG-USB2+ amplifier (OTBioelettronica s.r.l., Torino, Italy), using OTBiolab Software (OTBioelettronica s.r.l., Torino, Italy). Disposable Silver/Silver chloride surface electrodes with conductive gel (HEX Dual Electrodes, Noraxon USA, Scottsdale, AZ, United States) were placed on the skin of the subject. A moistened conductive band was wrapped around the wrist, ensuring contact with the radial and ulnar stylar processes, to serve as a reference electrode. Bipolar cables were attached to the disposable electrodes and connected to the amplifier.

Force Measurement
An ATI Mini 27 Ti force sensor integrated in the MR-StretchWrist is used to measure the wrist flexion/extension torque (Full-scale Load (FSL): 2 Nm, Resolution: 0.5 mNm, Measurement uncertainty: 1.5% of FSL). Transducer signals are processed by a preamplifier box (ATI Industrial Automation, Apex, NC, United States), and digitized by an analog/digital I/O device (PCIe-6321, National Instruments, Austin, TX, United States) connected to a laptop.

Control Software
Software for robot position control, perturbation timing, and for the graphical user interface was developed in Matlab and Simulink, and executed in real-time (sample rate: 1 kHz) using the QUARC real-time control software (Quanser, Markham, ON, Canada). Encoder data were acquired using the Q2 USB data acquisition board (DAQ).

EMG Electrode Positioning
To determine the location of the electrodes to measure activity of the FCR and ECU, manual palpations were performed on the right arm of the subjects. Repeated wrist flexion during palpation while the wrist was in its neutral position aided in locating the FCR. Similarly, repeated wrist extension in this same orientation aided in locating the ECU. Points parallel to the muscle's fibers were drawn 3 cm apart. The skin was then prepped with 70% isopropyl alcohol wipes and the application a thin layer of conductive skin prep gel (Nuprep, Weaver and Co., Aurora, CO, United States).

Experimental Procedures
In this study, subjects were seated with their forearms resting in a stationary support connected to a normal desk, with the forearm extending anteriorly in front of the body. Their hand was strapped inside of a mold such that their wrist was in a semiprone/neutral condition, and any static wrist flexion/extension torque (range: 0 -2 Nm) would be supported by the mold with little deflections. A GUI was shown on a computer screen indicating the amount of wrist flexion/extension torque to apply. After the subject reached the appropriate torque target, and maintained it within a 50 mNm range for 0.4 to 0.8 s, the robot perturbed the wrist in the direction that would stretch the agonist muscle with respect to the cued torque (i.e., if the background torque was flexion, the perturbation was wrist extension). Perturbations were applied for a duration of 200 ms in all conditions, so to avoid undesired oscillations of the EMG signal due to impact dynamics arising from the abrupt end of a perturbation. After each wrist perturbation, the robot halted for 1 second before returning to the neutral position for the following perturbation. Numerous conditions were studied in this protocol, defined by factorial combinations of four factors: (1) perturbation velocity, (2) perturbation direction, (3) background torque, (4) task instructions.
Perturbation velocity assumed three levels: 50, 125, and 200 deg/s. Perturbation direction assumed two levels: wrist flexion or wrist extension. Background torque was set to either 0 or 200 mNm. Task instruction assumed two levels: "yield" (Y) and "do not intervene" (DNI). In the Y condition, subjects were told to not provide any resistance after the perturbation and yield (i.e., relax) to the movement. In the DNI condition, subjects were told to continue applying the same amount of torque that they were applying prior to perturbation. These two instructions are fundamentally the same in the absence of a background torque.
Ten trials per condition (combination of each level of velocity, direction, background torque, instruction) were collected, resulting in each experiment consisting of a total of 240 perturbations. The order in which conditions were applied was randomly generated by the Simulink and MATLAB files. Furthermore, the time between the end of the prior perturbation and the start of the next perturbation cue was randomized between 3 and 7 s.

EMG Processing
Pre-processing of the raw EMG data was conducted using Matlab code (Matlab 2017a, MathWorks, Natick, MA, United States). After on-board analog amplification (gain = 2000, amplifier settings: low-pass cut-off frequency: 0.3 Hz, high-pass cut-off frequency: 4.4 kHz), data were sampled at 10240 Hz. A bandpass filter, a 4th order Butterworth filter (f LP = 250 Hz and f HP = 20 Hz), was used to remove low frequency noise related to movement artifacts and high frequency noise related to intrinsic measurement noise. Matlab's filtfilt function was used to obtain a filtered waveform without any phase shift of the signal. The filtered waveform was then rectified. Digital outputs produced by the DAQ board (Q2-USB, Quanser, Markham, ON, Canada) were sent to the EMG amplifier to identify instants of perturbation onset, used for segmentation of the EMG signal into 200 ms long timeseries, one for each perturbation, each starting at the time of perturbation onset. EMG tracks were shifted in time such that t = 0 ms corresponded to the instant of arrival of the pulse.
The amplitude of the segmented EMG signal measured during a perturbation was normalized by the average magnitude of the rectified EMG signal measured from that same muscle during agonist background contractions preceding all perturbations. This procedure allowed us to conduct group analysis of normalized EMG data. With this procedure, an EMG signal of unitary magnitude indicated that the rectified EMG signal measured during each perturbation had the same amplitude of the average EMG signal generated by that muscle for a 200 mNm isometric torque. Units of the normalized EMG timeseries are referred to as normalized units (nu).
These procedures resulted in 10 segmented timeseries extracted from each subject per combination of conditions (velocity, direction, torque, instruction). Long-latency response amplitude (LLRa), defined as the average signal of the EMG tracks in the [50 100] ms interval, was calculated for each perturbation for both muscles, and indexed as a function of subject, repetition, and combination of experimental conditions. The subject-specific average of LLRa for each combination of perturbation conditions was used as the outcome measure for the 0-D statistical analysis.
Repeated measurements from each subject were averaged to yield the timeseries with average rectified EMG response for each subject EMG sub,ν,d,t,inst used for the 1-D statistical analysis. Group averages EMG ν,d,t,i and corresponding standard deviations s ν,d,t,i were then calculated for display purposes.

Statistical Analysis
Two four-way full factorial linear mixed model ANOVAs were conducted using the subject-specific average LLRa measured for FCR and ECU, respectively, as outcome measure. The four factors included in the ANOVA were perturbation velocity (0, 125, 200 deg/s), background torque (0, 200 mNm), instruction (yield, do not intervene), and perturbation direction (stretch vs. shorten), defined based on the effect that the perturbation would have on the length of each muscle. As an example, a flexion perturbation would correspond to the "stretch" level for ECU, and the "shorten" level for FCR. Statistical analysis was conducted using JMP, and all variables were coded as nominal variables, with the default conditions of 50 deg/s for velocity, 0 mNm for background torque, yield for instruction, non-stretch for perturbation direction. JMP Pro Version 14 (SAS Institute Inc., Cary, NC, United States) was used for this analysis. The linear mixed model included 16 terms for fixed effects (4 main effects, 6 two-way interactions, 4 three-way interactions, one four-way interaction term, and an intercept), plus an additional set of offset variables for random subject-specific effects. The Satterthwaite method was used to determine the number of degrees of freedom in the model. All terms are reported if their estimated effect is significant at the type I error rate α = 0.05. In those cases, Tukey HSD post hoc tests were also used to determine pairs of levels with significant differences.
A 1-D ANOVA was also conducted on the timeseries of rectified EMG signal (EMG sub,ν,d,t,i ) measured during the postperturbation interval comprised between 0 and 200 ms using the spm1d software (Pataky et al., 2015(Pataky et al., , 2016. 1-D statistical analysis models are useful to analyze the effects of the experimental conditions on the perturbation-induced muscle response without prior hypotheses on the specific time interval where an effect is expected (Pataky et al., 2015(Pataky et al., , 2016. While with the 0D analysis we restricted our focus on the time interval ensuing the perturbation comprised between 50 and 100 ms (thus obtaining a scalar, or 0D, outcome measure), with the 1D model we sought to determine whether there is an effect on the timeseries of measured EMG amplitudes associated with all the experimental conditions, and their combinations, at any time point in the postperturbation interval comprised between 0 and 200 ms. As such, the 1-D ANOVA is useful to establish effects of all factors and their interaction at multiple time-points, but controlling for the multiple comparisons resulting from this type of analysis.
Because the current version of the spm1D Matlab software only allows to build full factorial models with a maximum of three main effects, we broke each four-way model into two three-way models and performed our analysis using the spm1d function anova3rm. Each model included the factors speed, background torque, and instruction, with one model including LLRs measured during stretch, and the other model including LLRs measured during shortening. This setup allowed us to study the effects and interaction of all factors studies in the 0D analysis, with the exception of all terms involving perturbation direction. Given the software limitations, we chose to exclude perturbation direction from this analysis as it was an effect that had been largely neglected by most other studies.
The 1-D analysis was implemented as a mixed model ANOVA for measurements collected at all time points, which controls for the associated multiple comparisons using random field theory (Pataky et al., 2016). As such, the output of the 1D ANOVA procedure is a time-series of F scores for all main effects and their interaction, combined with the identification of time intervals where those effects are significant at α = 0.05, corrected for the multiple comparisons performed at multiple time points.
When an effect or an interaction was significant during the LLR time-window (i.e., from 50 to 100 ms after perturbation onset), we conducted post hoc tests to determine pairs of levels with significant modulation in EMG signal.

RESULTS
Of the 13 participants recruited for this study, data sets of 2 subjects for ECU and 3 subjects for FCR were excluded from analyses due to poor EMG recordings. Therefore, statistical analysis was performed on data collected on n = 11 individuals for the ECU muscle, and n = 10 individuals for the FCR muscle. Timeseries extracted in the different experimental conditions for the ECU muscle are shown in Figures 2, 3. Similar representations are provided for the FCR muscle in Supplementary Figures 1, 2

0-D Analysis
The linear mixed model computed an adjusted R 2 of 0.627 for FCR, and an adjusted R 2 of 0.788 for ECU. The model reported a significant effect of all four main factors. However, since all factors are involved in several two-and three-way interactions, only the interactions will be analyzed and discussed below. Results are presented below as least square means ± standard error (LLRa) or difference in least square means ± standard error FIGURE 2 | Time series of EMG signal (left half-panels) and LLRa (boxplot in the right half-panels) from all subjects resulting with perturbations in the flexion direction (stretching the ECU muscle). Color indicates the speed of the perturbation (red = 50 deg/s, orange = 125 deg/s, green = 200 deg/s), with the line indicating the mean and the shaded area indicating the ± 1 s.d. region. Panels in different rows include measurements at two levels of task instruction -Yield (top), and DNI (bottom); columns include measurements at two levels of background torque -200 mNm (right), and 0 mNm (left). In the boxplots, darker shades represent 1 standard deviation from the mean, and lighter shades indicate the 95% confidence interval. ( LLRa) in units of the outcome measure, i.e., normalized units (nu). A report including the significant fixed effects from the linear mixed model ANOVA is provided below in Tables 1, 2. A significant three-way interaction between perturbation direction, instruction, and torque was measured for the ECU muscle (Figure 4). A significant effect of instruction or perturbation direction on LLRa was observed only in presence of a background torque for both muscles. In presence of background torque, the LLRa associated with stretch perturbations was greater than the LLRa associated with shortening perturbation both when task instructions were DNI (LLRa -Stretch DNI: 3.80 ± 0.17 nu, Shorten DNI: 0.79 ± 0.17 nu, p < 0.001), and when task instructions were to yield (LLRa -Stretch Y: 2.04 ± 0.17 nu, Shorten Y: 0.74 ± 0.17 nu, p < 0.001). In presence of background torque and muscle stretch, LLRa was greater in the DNI condition compared to yield ( LLRa -DNI vs. Y: 1.76 ± 0.16 nu, p < 0.001). A significant three-way interaction between perturbation direction, torque, and velocity was measured for the ECU muscle ( Figure 5). A significant effect of perturbation velocity or direction was observed only in presence of a background torque. In presence of background torque, the LLRa associated with stretch perturbations was greater Five two-way interaction terms were significant for ECU, while three terms were significant for FCR. Three terms were common to both muscles and are described next. The interaction between instruction and torque resulted from a greater increase in LLRa measured in the DNI conditions compared to the yield conditions measured in the presence of a background torque for both muscles (Figure 6). In the presence of background torque, LLRa measured during the DNI condition was significantly greater than LLRa measured in the Y condition (FCR LLRa -DNI: 3.67 ± 0.53 nu, Y: 2.39 ± 0.53 nu, p = 0.001; ECU LLRa -DNI: 2.30 ± 0.15 nu, Y: 1.40 ± 0.15 nu, p < 0.001). There was no difference in LLRa measured in absence of background torque between the two instructions (FCR LLRa -DNI: 1.12 ± 0.53 nu, Y: 1.22 ± 0.53 nu, p = 0.990; ECU LLRa -DNI: 0.52 ± 0.15 nu, Y: 0.52 ± 0.15 nu, p = 1). Significant differences between torque levels were measured in both the FCR and ECU for both instructions (FCR LLRa -DNI: 2.55 ± 0.34 nu, p < 0.001, Y: 1.16 ± 0.34 nu, p = 0.004; ECU LLRa -FIGURE 6 | Least square means for all significant two-way interactions for both the FCR (A) and the ECU (B). Letters on the plots are representative post hoc Tukey HSD tests; pairs of elements that do not share a letter are significantly different. Post hoc tests are, and the corresponding letter notation are separate for each panel.
The interaction between perturbation direction and velocity resulted from a greater increase in LLRa measured when muscles are stretched by a perturbation compared to when they are shortened measured at higher velocities (Figure 6) The significant interaction between perturbation direction and torque resulted from the greater increase in LLRa associated with stretch perturbations compared to shortening measured in presence of background torque (Figure 6). LLRa associated with stretch perturbations were greater than those associated with shortening in presence of background torque (FCR LLRa -0 mNm: 1.62 ± 0.34 nu, 200 mNm: 3.95 ± 0.34 nu, p < 0.001; LLRa ECU -0 mNm: 0.27 ± 0.11 nu, 200 mNm: 2.15 ± 0.11 nu, p < 0.001). LLRa measured in presence of background torque were greater than in absence of background torque for both muscles when they were stretched (FCR LLRa stretch: 3.02 ± 0.34 nu, p < 0.001; ECU LLRa stretch: 2.26 ± 0.11 nu, p < 0.001), but only for ECU when shortened (FCR LLRa shorten: 0.69 ± 0.34 nu, p = 0.180; ECU LLRa shorten: 0.38 ± 0.11 nu, p = 0.006).
The second two-way interaction term significant only for the ECU muscle was the interaction between perturbation direction and instruction. This term resulted from a greater increase in LLRa measured in the DNI condition compared to yield condition measured when muscles are stretched (Figure 6). When muscles were stretched, LLRa associated with the DNI conditions were larger than the yield condition (LLRa -DNI: 2.23 ± 0.15 nu, Y: 1.35 ± 0.15 nu, p = 0.001), while no significant difference between instruction conditions was measured when muscles were shortened. A significant increase in LLRa was measured in both instruction conditions when muscles were stretched compared to shortened, but this increase was greater in the DNI condition compared to the yield condition ( LLRa -DNI: 1.64 ± 0.11 nu; Y: 0.79 ± 0.11 nu, p < 0.001).

1D Analysis
The results of the 1D ANOVA are represented in terms of a timeseries of F scores, shown in Figures 7, 8 for the ECU and FCR, respectively. For effects and interactions that have significant upcrossings within the LLR region, post hoc 1D t-tests were used to further break down the effects.
The three-way interaction between instruction, torque, and velocity was not significant within the LLR for both muscles and both directions, however, there is a narrow significant upcrossing during the SLR for ECU Stretch (25.8 to 28.3 ms, peak F 2,20 = 11.495).
The two-way interaction between instruction and torque is significant within the LLR region for both directions of the ECU muscle, but not for the FCR muscle. Within the LLR region, there are significant upcrossings for ECU during both muscle shortening and stretching (ECU Stretch: 52.6 to 70.0 ms, local peak in LLR region F 1,10 = 26.152; ECU Shorten: 62.2 to 66.8 ms, peak F 1,10 = 37.434). There are also significant upcrossings within the voluntary region for the stretch of the ECU and FCR (ECU Stretch: 148.6 to 200 ms, peak F 1,10 = 26.152; FCR Stretch: 183.0 to 196.0 ms, 199.3 to 200 ms, peak F 1,10 = 27.377).
The two-way interaction between instruction and torque measured for the ECU muscle is broken down in the 1D post hoc t-tests to analyze the measured effect at all time points, as shown in Figure 9. Analysis of post hoc tests highlights how the DNI has positive or negative effect on the amplitude of processed EMG recordings at different time points, and the effect differs as a function of stretch and background torque condition. When the ECU was stretched, the normalized EMG signal in DNI was greater than yield within the LLR and voluntary regions in the presence of background torque, (ECU Stretch: 49.1 to 200 ms, peak T = 13.703, Figure 9B, center). No significant change in EMG signal was measured in absence of background torque as a function of task instructions (Figure 9B, left). As a result, the change in EMG signal measured between the DNI and Y condition was greater in presence of background torque than in absence only after the delay similar to that considered for forearm LLRs (ECU Stretch: 38.7 ms to 200 ms, peak T = 14.450, peak in the LLR region T = 9.397, Figure 9B, right). When the ECU was shortened, EMG recordings measured in presence of background torque in the DNI condition were greater than those measured in the yield condition in the initial part of LLR, but then EMG recordings measured during DNI were smaller than yield at a later time in the LLR time period (ECU Shorten: 57.7 to 77.0 ms, local peak T = 4.470, 94.6 ms to 129.5 ms, local peak T = -3.925, Figure 9A, center). Instead, no significant change in EMG signal was measured in absence of background torque as a function of task instructions (Figure 9A, left). The change in T scores between torque levels for the shortened condition also indicated significant upcrossings in the LLR and voluntary regions (ECU Shorten: 54.2 to 73.5 ms, local peak T = 4.961, 92.5 to 132.9 ms, local peak T = −5.461, Figure 9A, right).
The two-way interaction between torque and velocity is significant in the LLR region for the ECU when stretched and for the FCR when it is both shortened and stretched. For both the ECU and FCR there is a significant upcrossing within the LLR region (ECU Stretch: 52.6 to 74.7 ms, peak F 2,20 = 19.872; FCR Shorten: 47.9 to 65.8 ms, peak F 2,20 = 22.421; FCR Stretch: 58.0 to 68.6 ms, peak F 2,20 = 19.244).
The two-way interaction between torque and velocity for muscle and perturbation direction is broken down in 1D post hoc tests for ECU Stretch, FCR Shorten, and FCR Stretch (Figures 10-12). For the stretch of the ECU (Figure 10), significant upcrossings were present in post hoc tests comparing EMG signals measured at different velocity levels and at multiple torque levels, primarily during the LLR time period. Analysis of the timeseries of t-scores demonstrate that EMG signal increases with velocity and with background torque, and that the region where such an increase is measured largely overlaps with the expected latency of a LLR. The largest upcrossing was measured for the 200-50 deg/s comparison in presence of a background torque (40.8 to 98.1 ms, peak T = 10.202, Figure 10, center row, center column). Upcrossings within the LLR region were also measured for the 200-50 deg/s t-test in the absence of background torque, the 200-125 deg/s t-test in the presence of background torque, and the 125-50 degree/second t-test in both torque conditions. Significant upcrossings were measured within the LLR region for the 200-50 deg/s comparison and the 200-125 degrees/second comparison ( 200-50 deg/s: 56.4 to 73.8 ms, peak T = 5.915, Figure 10 center row, right column, 200-125 deg/s: 56.0 to 62.7 ms, peak T = 4.304, Figure 10, top row, right column). The significance of the between-torque condition difference of the t-scores resulting from comparing pairs of velocity levels indicates that background torque differentially modulates the velocity dependence of EMG signals.
Qualitatively similar results were measured for the stretch of FCR (Figure 11). Significant upcrossings were measured for t-test comparisons mostly during the LLR time period. The largest t-scores were generated for the 200-50 deg/s t-test in the presence of background torque (23.7 to 80.9 ms, peak T = 12.856, Figure 11 center row, center column). Significant upcrossings were also present for the 200-125 deg/s t-test in presence of background torque, 200-50 degree/second t-test in the absence of background torque, the 200-125 deg/s t-test in both torque conditions, and the 125-50 deg/s t-test in both torque conditions. One significant upcrossing was measured for the between-torque condition difference of t-scores resulting from comparing pairs of velocity conditions ( 200-50 deg/s: 57.9 to 64.0 ms, peak T = 4.604, Figure 11 center row, right column). For FCR responses measured during muscle shortening (Figure 12), significant upcrossings were present in t-test comparisons between velocity levels in the presence of background torque. Like in the previous conditions, the largest t-scores were generated for the 200-50 degree/second t-test in the presence of background torque (46.0 to 93.3 ms, peak T = 8.601, Figure 12 center row, center column). There are also smaller, significant upcrossing for this t-test in the voluntary region (149.1 to 165.1 ms, 187.9 to 200 ms). Significant upcrossings were also present for the 200-125 deg/s t-test in presence of background torque, 200-50 degree/second t-test in the absence of background torque, and the 125-50 degree/second t-test in both torque conditions. Significant upcrossings within the LLR region were calculated for difference between t-tests of different background torque levels ( 200-50 deg/s: 48.7 to 63.2 ms, peak T = 6.587, Figure 12 center row, right column, 200-125 deg/s: 50.6 to 59.2 ms, peak T = 5.176, Figure 12, top row, right column). Such differential effect of background torque at multiple velocities on EMG signal collected from a muscle under shortening was only observed for the ECU muscle in the 0D analysis (Figure 5 The two-way interaction between instruction and velocity is only significant for the shortened state of the FCR during a narrow time window outside of the LLR region (189.8 to 193.1 ms, peak F 2,20 = 10.959).
With respect to the main effects, both the FCR and ECU had significant upcrossings within the LLR region. The main effect of instruction was significant within the LLR region for the stretch of the FCR and ECU only for a very short time period (FCR stretch: 69.8 to 73.2 ms, local peak F 1,10 = 26.080; ECU stretch: 81.1 to 82.8 ms, local peak F 1,10 = 22.546) and has upcrossings within the voluntary region for both the stretch of the FCR and ECU (FCR Stretch: 184.6 to 200 ms, peak F 1,10 = 30.620; ECU Stretch: 151.1 to 200 ms, peak F 1,10 = 32.771). There were otherwise no significant upcrossings for the effect of instruction for the shortened state of the muscles.

DISCUSSION
In the present study, we aimed to investigate the effects of and the interactions among four behavioral factors in modulating the LLR amplitude (LLRa) evoked from FCR and FIGURE 9 | Results of the 1D t-tests for the interaction between task instruction and torque for the shortened (A) and stretched (B) states of the ECU. Plots indicate the statistical score of 1D paired t-tests between the DNI and yield instructions and split in columns by torque level. Plots in the right column report the difference between the T-scores in the second and first column's tests, used to demonstrate interaction between these factors. The dashed red line indicates the Bonferroni-adjusted critical T values and regions in green indicate significant upcrossings.
FIGURE 10 | Results of the 1D t-tests of the interaction between torque and velocity for the stretched state of the ECU. Plots in different rows indicate the statistical score of 1d paired t-tests between the velocity levels (200-125, 200-50, 125-50 deg/s) and columns indicate different torque levels. Plots in the right column report the difference between the T-scores in the second and first column's tests, used to demonstrate interaction between these factors. The dashed red line indicates the Bonferroni-adjusted critical T-values and regions in green indicate significant upcrossings.
ECU muscles during ramp-and-hold perturbations applied to the wrist joint. The four behavioral factors studied in this work were perturbation direction (stretch vs. shorten), background muscle activation (0 vs. 200 mNm of joint torque requiring agonist muscle activity), perturbation velocity (50, 125, and 200 deg/s), and task instructions (yield vs DNI). In line with previous studies (Lee and Tatton, 1982;Bedingham and Tatton, 1984;Calancie and Bawa, 1985;Miscio et al., 2001;FIGURE 11 | Results of the 1D t-tests of the torque and velocity interaction for the stretched state of the FCR. FIGURE 12 | Results of the 1D t-tests of the torque and velocity interaction for the shortened state of the FCR. Lewis et al., 2004Lewis et al., , 2005Lewis et al., , 2006Lewis et al., , 2010Schuurmans et al., 2009;Kurtzer et al., 2014), we observed that all factors modulate the LLRa in the FCR and ECU muscles. However, given the use of full factorial design involving all combinations of the four factors, we were able to study higher-order interactions among the four behavioral factors. Moreover, the use of a 1D analysis of stretch-evoked EMG data allowed us to establish the effects of the behavioral factors on EMG amplitude at multiple timepoints in the post-perturbation period, without restricting the analysis to a-priori selected time windows where a modulation would be expected.
In line with previous studies (Tatton and Bawa, 1979;Bedingham and Tatton, 1984), our findings demonstrated that perturbation velocity significantly modulates LLRa, as larger LLRa resulted from perturbations with 200 deg/s velocity compared to the 125 and 50 deg/s velocity conditions, while the 50 deg/s perturbation velocity evoked the lowest LLRa. Note that given the chosen task design, where perturbations are applied for a constant duration (200 ms), perturbation velocity and displacement are associated, so the measured effect could be due to the fact that the joint is subject to a greater displacement in the higher velocity conditions. From studies that systematically spanned a set of perturbation velocity and durations (Lee and Tatton, 1982;Lewis et al., 2005;Schuurmans et al., 2009), we know that total joint displacement is the primary contributor to LLRs, with perturbations at high velocity but brief duration not generating significant LLRa. In our experience, the execution of LLRa studies on perturbation durations that require stopping the limb within the [50 100] ms interval proved to be challenging, as in preliminary data we observed the occurrence of motion artifacts associated with the oscillations induced by the manipulator stopping. For this reason, we decided to focus on perturbations of longer duration (200 ms), which allowed to completely avoid those artifacts.
Results also revealed a significant increase in LLRa when a 200 mNm background torque was applied prior to the perturbations in comparison to the perturbations with no background torque. Pre-existing background muscle activation is thought to reflect an automatic adjustment mechanism, known as the automatic gain component of LLR (Bedingham and Tatton, 1984;Matthews, 1986;Miscio et al., 2001;Pruszynski et al., 2009).
Task instruction also significantly modulated the LLR in this study. The majority of previous studies examined the 'yield' or 'DNI' instruction with a 'resist' or 'compensate' instruction (Miscio et al., 2001;Lewis et al., 2006;Kurtzer et al., 2014). However, no prior study tested whether the 'yield' and 'DNI' task instruction differentially modulate the LLRa, with the exception of one study by Calancie and Bawa (1985) who recruited only two healthy individuals. Our findings showed that larger LLRs were evoked when participants attempt keep constant the muscles' activation state (i.e., DNI), compared to when they are asked to yield to the perturbations, which is justifiable in accordance to the evidence that the temporal overlap of two different responses including a task-dependent response and an automatic response results in the task-dependent change in LLR amplitude (Rothwell et al., 1980;Lewis et al., 2006;Pruszynski et al., 2011).
Also, in line with a few previous studies quantifying LLR from both stretched and shortened muscles (Miscio et al., 2001;Lewis et al., 2004), the above-mentioned modulations had more prominent effects on the LLRa evoked from the muscle stretched in response to the perturbations vs. the shortened muscle, however several significant effects of the behavioral factors were also measured in muscles undergoing shortening (see below for details).
With our full factorial design, we were able to observe more granular modulations of the LLR response induced by the four behavioral factors reported above. Specifically, significant results were observed for all the three-way interaction terms that can be computed without including the interaction between instruction and velocity, as detailed below.

Perturbation Direction, Task Instruction and Background Torque
As illustrated in Figure 4, perturbation direction did not affect the LLRa evoked from the perturbations with no background torque (i.e., the LLR measured during perturbations shortening the muscle were of similar amplitude as the one measured during perturbations stretching the muscle). This finding is in contrast with a previous study by Miscio et al. (2001), which demonstrated that there would be a lower LLRa evoked from the shortened muscle in comparison to the LLRa evoked from the stretched muscle even in absence of background muscle activity.
In addition, task instruction did not modulate the LLRa evoked in absence of background torque. This result is expected and serves as a validation of the pursued experimental design, as in this case the behavior elicited in the Y and DNI conditions is expected to be identical.
Results from the perturbations with 200 mNm background torque showed a significant difference between LLRa evoked from stretched muscle in the Y and DNI task instructions. Larger LLRa was detected from perturbations when participants were instructed to maintain their level of activation (DNI) after the perturbation, compared to when they were instructed to yield to the displacement. In contrast, task instructions did not modulate the LLRa evoked from the shortened muscle neither in presence of a 200mNm background torque. Also, the LLRa for shortened muscle was significantly lower than the LLRa for stretched muscle in both instruction conditions. This interaction is also consistently reported also in the 1D analysis (Figure 9).
In agreement with this result, Miscio et al. (2001) also reported a smaller LLRa for shortened muscle in comparison to the stretched muscles when there was a 10% maximum voluntary contraction (MVC) background torque and subjects were instructed to not intervene to the perturbations. However, their task instructions did not include a 'yield' condition for comparison. The EMG activity recorded from the shortened muscle in the previous study was interpreted as a volume conducted response based on their direct nerve stimulation results. Besides, their regression analysis showed a positive correlation between the area of FCR stretch response and the low-level activity recorded from the shortened antagonist muscle (ECR) (r 2 = 0.59). Instead, our data supports a differential modulation of the response for the stretched and shortened muscles between different levels of task instructions, suggesting that the source of the observations may not be due to cross-talk, but to true physiological decoupling of the muscle responses.

Perturbation Direction, Background Torque, Perturbation Velocity
Results of the 0D analysis ( Figure 5) support a significant interaction between the effects of torque and velocity on the LLRa evoked from perturbations stretching the ECU muscle, but not from perturbations shortening the ECU muscle. In support of this observation, the 1D analysis (Figure 10) revealed a significant interaction between torque and velocity on the LLRa for ECU muscle only under stretching, but not under shortening. For the FCR muscle, the interaction was significant in response to the both stretching and shortening perturbations, though the differential effects of background activity are smaller, and last for a smaller duration for this muscle (Figures 11, 12). The smaller size and duration of the effects of torque and velocity, combined with the greater between-subject variability of FCR data might be a possible reason why the 0D analysis did not detect a significant effect of this interaction term for the FCR muscle.
Previous studies also demonstrated a significant effect of the interaction between perturbation torque and velocity on the LLRa evoked from the stretched muscle (Bedingham and Tatton, 1984). However, based on our literature review, no prior study evaluated this interaction for the muscle subject to a perturbation that shortened the muscle. Berardelli et al. (Berardelli and Hallett, 1984), reported that background torque and velocity would increase the amplitude of the EMG activity evoked from the shortened tibialis anterior muscle with a latency within 100-150 ms, but no interaction analysis was conducted in that study.
For both muscles, perturbation direction did not modulate the LLRa elicited by any velocity condition in absence of background torque. However, Miscio et al. (2001), observed lower LLRa for the shortened muscle compared with the stretched muscle which were perturbed with 500 deg/s velocity with no background torque. It might be possible that higher perturbation velocities are needed to differentiate the long latency stretch response of the agonist and antagonist muscles when there is no background activity. Moreover, to the best of our knowledge, our study is the first study comparing the effects of different perturbation velocities on the LLRa between the shortened and stretched upper limbs' muscles, and the first to observe an interaction between background torque and velocity in the LLRa of a muscle undergoing shortening, an effect that was visible for the FCR muscle.

Interaction Between Task Instruction and Velocity
No analysis has revealed any interaction between the task instruction and perturbation velocity. The majority of previous work studied the interaction between perturbation velocity with perturbation duration and amplitude (Lee and Tatton, 1982;Lewis et al., 2005Lewis et al., , 2006Schuurmans et al., 2009). Among them, Lewis et al. (Lewis et al., 2006), was the only study that examined the interaction between task instruction and perturbation velocity. In contrast to our results, they reported a significant interaction, which resulted in a greater facilitation in the LLRa evoked from the Biceps brachii muscle, for the 'Flex' instructed perturbations with a 90 degree/s velocity. In their 'Flex' instruction, subjects were instructed to flex their elbows as soon as possible in response to the perturbation, a condition which wasn't included in our experiment. Hence, it is possible that countering perturbations applied to the proximal upper limb joints (such as the elbow) would elicit a different modulation in the LLR from the one observed here for forearm muscles.

Similarity of the Effects Between FCR and ECU
Some of the interaction terms were only significant for ECU, and not for FCR. Specifically, no three-way interaction term was significant for the FCR muscle (as opposed to two terms significant for FCR), and two of the five two-way interactions significant for the ECU muscle were not significant for FCR. Overall, the mismatch of results between the two muscles is primarily due to larger variability of the measurements obtained for the FCR muscle compared to the ECU muscle, while all effects were generally in the same direction for both muscles. In many cases, in fact, the magnitude of measured modulations in LLRa (expressed in normalized units, so roughly comparable between muscles) was actually larger for the FCR muscle than for the ECU muscle. However, the standard error of the mean estimated by the mixed model was several times larger for the FCR muscle compared to the ECU muscle.
As an example, the term corresponding to the three-way interaction between direction, instruction, and torque was close to the conventional type-I error rate for the FCR muscle (p = 0.055). Further inspection of the parameter estimates corresponding to that interaction term shows that the mean difference in LLRa measured between Y and DNI conditions was slightly larger in the FCR muscle than it was for the ECU muscle under stretch (FCR LLRa -DNI: 6.05 ± 0.58 nu, Y: 3.95 ± 0.58 nu; ECU LLRa -DNI: 3.80 ± 0.17 nu, Y: 2.04 ± 0.17 nu), while smaller effects were measured in the shorten condition (FCR LLRa -DNI: 1.28 ± 0.58 nu, Y: 0.83 ± 0.58 nu; ECU LLRa -DNI: 0.78 ± 0.17 nu, Y: 0.73 ± 0.17 nu) for both muscles. Yet, the four-fold larger s.e.m. measured for the FCR prevented this effect from becoming significant.

0D vs. 1D Analysis
To the best of our knowledge, this study is the first to apply a mixed-model 1D analysis to the timeseries of EMG signals measured in response to velocity-controlled perturbations. 1D analyses are used to quantify if and when a timeseries is significantly modulated by an experimental factor, without a-priori hypotheses on the duration of this modulation, or on the time window where such a modulation is expected. Typically, reflex studies use 0D variables such as the average of the rectified signal in a pre-defined time window (Kurtzer et al., 2008), or the cumulative sum of the rectified EMG signal within a region of supra-threshold response (Brinkworth and Türker, 2003) -the latter definition allowing to explicitly account for a combined measure of amplitude and duration of a reflex response.
A shortcoming of using 0D datapoints generated from a 1D data set is the possibility for false positives. It has been indicated that the false positive rate of the 0D data is greater than the desired false positive rate of α = 0.05 when using noisy outcome measures such as peak in an interval (Pataky et al., 2016), and that this may be the case also in presence of smoothing (Pataky et al., 2015(Pataky et al., , 2016. Less noisy outcome measures such as the mean EMG signal in a pre-defined time window are not likely to result in limited control of the false positive rate of inference tests. However, the main advantage of using a 1D analysis for the study of stretch reflexes is the fact that the 1D analysis gives a more granular insight on the shape of the waveform (e.g., one peak, two peaks) which can be masked by the operation of averaging necessary for the creation of a 0D dataset. Moreover, the 1D analysis allows to test hypothesis on whether an effect is only present in a specific time-region, while controlling for the type-I error rate of inferential statements.
As an example, our 0D analysis showed that, in presence of background activation, there is an effect of task instruction on the LLRa measured from the ECU when stretched, but no effect measured when the muscle is shortened (Figure 4). However, a closer inspection on this modulation is provided by the post hoc t-tests presented in Figure 9. It is true that the difference between rectified EMG signal measured in the DNI and Y condition for the ECU under stretch is greater in the 200 mNm condition, reflecting a larger signal measured in the DNI condition (Figure 9, center). However, the same contrast plot extracted from the ECU during shortening perturbations highlights the presence of a biphasic response, where the DNI condition results in an early increase in EMG signal in the LLR time window, followed by a decrease in EMG signal. The average of the positive and negative changes is likely the reason why the same effect could not be captured using the 0D analysis.
Also, insight from a 1D analysis provides information about the gradient in timing resulting from different behavioral factors. Specifically, the effect of instruction is usually visible later than the one of velocity and torque, which have upcrossings in SLR and early LLR (Figures 7, 8). Instruction has a narrow region of significant main effects in the late LLR period (70 to 73 ms for FCR stretch, 81 to 83 ms for ECU stretch), with significant upcrossings in the voluntary period. Velocity conditions are significant in the early to mid LLR period (FCR stretch ∼40 to 70 ms, ECU stretch ∼45 to 90 ms, FCR shorten ∼50 to 90 ms, ECU shorten ∼45 to 80 ms), all having significant peaks near 60 ms. The 1D analysis allowed us to deduct the presence of a 10 to 20 ms delay between the significant peaks due to instruction versus velocity. The time frame of these significant peaks are comparable with previous findings and the identification of two distinct peaks within the LLR period, named M2 and M3 in the literature. However, it is difficult to compare the exact timing in absolute terms since "latency" in the literature is often defined relative to the onset of EMG activity (Lee and Tatton, 1982;Calancie and Bawa, 1985;Miscio et al., 2001;Lewis et al., 2006;Kurtzer et al., 2014).
Yet, the results of the 0D and 1D analyses do not perfectly overlap. This is in part also due to limitations in the current version of the software used for the analysis. The spm1D MATLAB package allows up to three input conditions, as seen through its anova3rm function. To analyze the four conditions studied, two analyses per muscle were conducted when in the shortened and stretched conditions. As such, none of the threeway interactions that were significant in the 0D analysis could be studied in the same form using the 1D analysis. Another limitation of the 1D analysis is related to its limited statistical power. Because of the large number of datapoints included in the analysis, and the fact that the EMG signal is somewhat noisy and composed of multiple resolution elements in a time window of 200 ms following a perturbation, 1D analyses afford a reduced statistical power compared to a 0D analysis that only looks at the average within a pre-defined time window (Pataky et al., 2015).

Study Limitations
The sample size available for this study is in line with the one used in previous behavioral experiments with LLR of forearm muscles (Miscio et al., 2001;Pruszynski et al., 2009;Lewis et al., 2010) in terms of number of subjects recruited. The number of datapoints collected from each subject is also in line with most previous studies, with 240 datapoints per subject collected. However, given the need to estimate a larger number of model terms, a larger sample size would have allowed more precise estimation of all model effects and reduced the variance of the outcomes between the different muscles and the different statistical analysis techniques pursued.
The estimate of the effects of background torque is potentially confounded by the fact that we only applied perturbations in one direction, and the direction was predictable only in the 200 mNm condition. Specifically, in the 200 mNm condition, we applied perturbations that would stretch the agonist muscle, and as such the direction was predictable by the subject. Instead, in the 0 mNm condition, perturbations were applied in a random direction, not predictable by the subject. As such, the coefficients estimated for the two levels of the factor background torque might be different also simply because in the 200 mNm condition, subjects were able to predict perturbation direction, while in the rest condition, they were unable to do so. It is possible that prior knowledge of perturbation direction modulates the LLRa. As such, it is possible that our model estimated the compound effect of two distinct factors: a factor due to muscle pre-load, and another factor due to knowledge of the perturbation direction. This ambiguity could be addressed in future studies by applying perturbations that either stretch or shorten agonist muscles after pre-load.

CONCLUSION
In summary, this study has quantified the effect of four behavioral factors -background torque, task instruction, perturbation velocity and direction -on the long-latency response (LLR) amplitude evoked from the FCR and ECU muscles during rampand-hold perturbations applied to the wrist joint in the flexion and extension direction. Our analysis demonstrated that all of those factors modulate LLRa, and that their combination nonlinearly contribute to modulating the LLRa. Specifically, all the three-way interaction terms that can be computed without including the interaction between instruction and velocity significantly modulated the LLR. The interaction analysis suggested that higher background torque augmented the LLRa evoked from the stretched muscle when subjects are asked to maintain their muscle activation in response to the perturbations.
Besides, higher perturbation velocity increased the LLRa evoked from the stretched muscle in presence of a background torque. Also, our analysis identified significant modulations of LLRa also in muscles shortened by the perturbation, including an interaction between torque and velocity, and an effect of both torque and velocity. While a lot of the behavioral factors listed above nonlinearly contribute to modulating LLRa, we observed that the effects of task instruction and velocity do not combine more than linearly to modulate the LLRa.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article have been uploaded on a zenodo repository, https://zenodo.org/record/ 4641292#.YGEzVC2cZTZ.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The University of Delaware Institutional Review Board, protocol no. 1097082-6. The patients/participants provided their written informed consent to participate in this study.