Disentangling acceleration-, velocity-, and duration-dependency of the short-and medium-latency stretch re ﬂ exes in the ankle plantar ﬂ exors

Motorized assessment of the stretch re ﬂ ex is instrumental to gain understanding of the stretch re ﬂ ex, its physiological origin and to differentiate effects of neurological disorders, like spasticity. Both short-latency (M1) and medium-latency (M2) stretch re ﬂ exes have been reported to depend on the velocity and acceleration of an applied ramp-and-hold perturbation. In the upper limb, M2 has also been reported to depend on stretch duration. However, wrong conclusions might have been drawn in previous studies as the interdependence of perturbation parameters (amplitude, duration, velocity, and acceleration) possibly created uncontrolled, confounding effects. We disentangled the duration-, velocity-, and acceleration-dependence and their interactions of the M1 and M2 stretch re ﬂ ex in the ankle plantar ﬂ exors. To disentangle the parameter interdependence, 49 unique ramp-and-hold joint perturbations elicited re ﬂ exes in 10 healthy volunteers during a torque control task. Linear mixed model analysis showed that M1 depended on acceleration, not velocity or duration, whereas M2 depended on acceleration, velocity, and duration. Simulations of the muscle spindle Ia afferents coupled to a motoneuron pool corroborated these experimental ﬁ ndings. In addition, this simulation model did show a nonlinear M1 velocity- and duration-dependence for perturbation parameters outside the experimental scope. In conclusion, motorized assessment of the stretch re ﬂ ex or spasticity using ramp-and-hold perturbations should be systematically executed and reported. Our systematic motorized and simulation assessments showed that M1 and M2 depend on acceleration, velocity, and duration of the applied perturbation. The simulation model suggested that these dependencies emerge from: muscle-tendon unit and muscle cross-bridge dynamics, Ia sensitivity to force and yank, and motoneuron synchronization. NEW & NOTEWORTHY Previous research and de ﬁ nitions of the stretch re ﬂ ex and spasticity have focused on velocity-dependence. We showed that perturbation acceleration, velocity, and duration all shape the M1 and M2 response, often via nonlinear or interacting dependencies. Consequently, systematic execution and reporting of stretch re ﬂ ex and spasticity studies, avoiding uncontrolled parameter interdependence, is essential for proper understanding of the re ﬂ ex neuro-physiology.


Introduction
Reflexes are an important mechanism within human movement control to cope with external perturbations during daily living. Specifically, the stretch reflex is the rapid motor response to counteract an unexpected lengthening of a muscle. Unfortunately, an exaggerated stretch reflex, that is, hyperreflexia or spasticity, is often present in people with brain or neural injuries, such as cerebral palsy or spinal cord injury (1,2). Hyperreflexia contributes to the movement disorder observed in these people, limiting their functional independence.
Motorized assessment of the stretch reflex involves imposing a joint movement and measuring the subsequent response in muscle activity. The advantage of motorized above manual assessment is that the stretch perturbations can be precisely controlled and standardized (3). After a sudden muscle stretch, three consecutive responses can be observed in the electromyography (EMG) in the lower limb: the short-latency (M1), medium-latency (M2), and long-latency (M3) response (4). Motorized assessment is important to gain understanding of the stretch reflex, its physiological origin and to differentiate effects of neurological disorders.
The interdependence of the amplitude, duration, velocity, and acceleration parameters is important to consider when investigating the effect of perturbation characteristics. Regarding these four parameters, perturbation signals can only be designed based on three independent parameters. For example, Dietz et al. (11) investigated the velocity-dependence, but scaling of perturbation velocity was achieved by scaling acceleration, creating a potential confounder. All other studies that investigated the velocity-dependence also potentially had acceleration as confounder, as none reported the acceleration profile used (8)(9)(10)(12)(13)(14)(15). Similarly, studies on muscle spindle firing dynamics are subject to the same interdependence. For example, Blum et al. (19) observed a relationship of the Ia afferent response's dynamic index with stretch velocity and initial burst with stretch acceleration. However, the simulated stretch velocity and acceleration were varied with a perfect correlation, thus the observed relations cannot conclusively be linked to either velocity or acceleration. In general, wrong conclusions might have been drawn in previous studies regarding the amplitude-, duration-, velocity-, and acceleration-dependence of muscle spindle dynamics and subsequent M1 and M2 response.
The goal of this paper is to disentangle the duration-, velocity-, and acceleration-dependence and their interactions of the M1 and M2 stretch reflex in the ankle plantarflexors. Ramp-and-hold perturbations with the amplitude parameter as dependent variable are used to investigate M1 and M2. Therefore, the amplitude dependency is not investigated. Based on the previously reported dependencies, we hypothesize that the M1 response depends on velocity and acceleration. Moreover, we hypothesize that the M2 response depends on duration, velocity, and acceleration. The M2 duration-dependence has only been reported in the upper limb (17). To disentangle the perturbation parameters under investigation, 49 unique perturbation profiles are used. In addition, a biophysical simulation model of the muscle spindle Ia afferents (19) coupled to a motoneuron pool (20) was implemented. This simulation model was used to corroborate the experimental findings, investigate stretch reflex dependencies across an extended set of perturbation parameters and gain a physiological understanding of the observed dependencies. The outcome of this study aims to help understanding of the stretch reflex and to stress the importance of a sound perturbation profile design in future stretch reflex studies.

MATERIALS AND METHODS Participants
Ten volunteers with no history of neuromuscular disorders participated in the study: age 26.4 yr (SD 1.9), two women and eight men. The EEMCS ethics committee of the University of Twente approved the study under reference number RP 2018-58. All participants provided written informed consent.

Apparatus
Participants were seated on an adjustable chair with the right foot connected to a robotic manipulator fixed onto the chair frame, see Fig. 1. The foot connection to the manipulator used a rigid footplate and Velcro straps. The posture was controlled for by supporting the upper body and leg using the adjustable chair. The chair ensured that knee and hip angles were fixed at 150 and 120 , respectively. Both knee and hip were defined at 180 for a perfectly straight posture. The starting manipulator position was set at a 90 ankle angle, defined as the angle between shank and foot determined using a goniometer. The ankle and manipulator axes of rotation were visually aligned at the start of the experiment, minimizing knee translation due to the applied ankle rotations.
A one degree-of-freedom (DOF) manipulator (Moog, Nieuw-Vennep, The Netherlands) applied ramp-and-hold perturbations stretching the ankle plantarflexors in the sagittal plane. Ankle angle and velocity (i.e., angular velocity) were represented by the angular position and velocity of the footplate measured using the actuator's encoder. Ankle torque was measured with a torque sensor placed between the actuator and footplate. Angle, velocity, and torque were recorded at 2,048 Hz, all defined positive in dorsiflexion e u q r o T ) m N ( -3 Figure 1. Overview experimental setup. Participants were seated on an adjustable chair. A manipulator applied dorsiflexion, ramp-and-hold perturbations around the right ankle joint, while measuring the response in muscle activity. Participants were instructed to keep a constant background torque using a feedback screen. The feedback screen showed a (red) plantarflexion torque target around À3 ± 0.1 Nm and a (blue) smoothed history of the torque exerted by the participant.

Experiment Protocol
Participants were instructed to keep background torque constant throughout the experiment using a feedback screen, see Fig. 1. The feedback screen provided biofeedback of a 6 s history of the smoothed (moving average, 200-ms window) measured torque and a À3 ± 0.1 Nm torque target, that is, the participant exerted a plantarflexion torque. The torque task was used to ensure a constant background muscle activation, uncorrelated to changes in perturbation parameters. As task instruction can influence the stretch reflex response, participants were instructed to not respond to the perturbations, similar to Finley et al. (17). Moreover, participants were instructed to generate background torque as if rotating the ankle without using the upper leg. To support these instructions, the influence of the ramp-and-hold perturbations on the torque feedback was attenuated. A constant value was shown during each perturbation, equal to the torque value shown just before perturbation onset.
Stretch reflexes were elicited using 49 unique perturbation profiles. These 49 perturbation profiles were the combination of 2 acceleration levels (140 and 175 rad/s 2 ), 3 velocity levels (2.0, 2.5, and 3.0 rad/s), and 10 duration levels (30-75 ms with 5 ms steps), see Fig. 2. As a result, stretch amplitudes ranged from 0.031 rad (1.76 ) to 0.165 rad (9.45 ). The experimental scope was limited to avoid excessive muscle fatigue and participant loss of attention. As existing data sets already focused on acceleration and velocity (12,13,16,17), we opted to include a broad range of duration levels. Combining all levels would give 60 unique perturbations, however 11 of these perturbations had an infeasible combination of parameters. Specifically, for short duration stretches, high velocities cannot be reached given the chosen acceleration levels.
Perturbation profiles were designed to disentangle the duration, velocity, and acceleration parameters, see Fig. 2. The acceleration levels (140 and 175 rad/s 2 ) were taken directly from Finley et al. (17), as a clear M2 response was present in the ankle plantarflexors at these levels. The acceleration profile consisted of four smooth transition shapes (6 samples, sinusoidal) with a constant level in-between. The profile was scaled linearly to achieve the chosen acceleration levels. Moreover, changing the length of constant acceleration periods allowed to set velocity levels. The velocity levels (2.0, 2.5, and 3.0 rad/s) were chosen in a similar range to previous experiments reporting velocity-dependency of the ankle stretch reflex (10,12,13). Changing the length of constant velocity periods allowed to set duration levels. The range of duration levels (30-75 ms) was defined based on the durationdependency shown in the wrist (14,15). A small resolution of 5 ms was chosen as the duration effect for M2 has not been explored in the ankle before and nonlinear effects may exist (14,15).
The experiment consisted of 12 blocks with all 49 perturbation profiles elicited exactly once per block. The order of the perturbations was randomized for each block. The For each profile the duration, maximum velocity, and maximum acceleration parameters were fixed with the amplitude as dependent parameter. Detailed time series are shown for the two highlighted perturbation profiles (plus and circle). Right: commanded and measured angle, velocity, and acceleration for the two highlighted perturbation profiles. The corresponding ensemble-averaged stretch reflex EMG response of a single representative participant are also shown. The maximum acceleration parameter was set by scaling the (blue shaded) sinusoidal shape transitions. The maximum velocity parameter was set by elongating the (red shaded) period with maximum acceleration. The perturbation duration was set by elongating the (green shaded) period with maximum velocity. The stretch reflex EMG response shows the (yellow shaded) 10-ms M1 window and 20-ms M2 window used to quantify reflex activity. DF, dorsiflexion; PF, plantarflexion.
stretch reflexes were elicited with a random 3-to 5-s interval between each perturbation. The blocks were executed with a 2-min break between each block and a larger 5-min break between blocks 4-5 and blocks 8-9 to prevent fatigue.

Data Analysis
The correct execution of the torque task was checked for each stretch reflex during data analysis. This check was necessary as stretch reflexes were applied continuously, even when participants did not maintain the desired torque. The background torque was computed as the average torque over the 200-ms period before perturbation onset. All stretches with a background torque deviating more than ±0.2 Nm from the À3 Nm instructed were rejected from further analysis.
EMG signals were high-pass filtered (2nd-order, 5 Hz, Butterworth) and rectified. For each muscle of every participant, the M1 and M2 analysis windows were determined via visual inspection. The M1 and M2 windows were set using the ensemble average of all stretch reflex responses of a participant, aligned at perturbation onset. The M1 analysis window was set centered around peak M1 activity with a 10-ms window width. A narrow window width was used as M1 timing was quite consistent across reflexes and to avoid contamination with M2 activity (17). Across subjects, the SOL and GM/GL M1 windows were placed starting at 39-53 ms (49-ms median) and 39-49 ms (47-ms median) after perturbation onset, respectively. The M2 analysis window was centered around peak M2 activity with a 20-ms window width. Contrary to M1, a wider 20-ms window was used for M2 to reflect the larger variability in timing observed compared with M1. The SOL and GM/GL M2 windows started at 54-70 ms (64-ms median) and 52-67 ms (62-ms median) after perturbation onset, respectively. Figure 2 depicts this difference between M1 and M2 timing for a representative participant.
For each stretch reflex of every participant, background EMG activity as well as M1 and M2 response magnitude were quantified. Background EMG should reflect an average activity over the period before perturbation onset. As a result, the background EMG was computed as the mean EMG activity over the 100-ms period before perturbation onset. M1 and M2 response measures should reflect the true reflexive magnitude, typically appearing as a double-peak shape due to rectification. To compensate for background activity, background EMG was subtracted from the reflexive response and the resulting signal was half-wave rectified. M1 and M2 magnitudes were defined as the root mean square (RMS) value of this half-wave rectified signal within the M1 and M2 analysis windows, respectively. Finally, for each perturbation profile the background EMG, M1 magnitude, and M2 magnitude measures were averaged across all repetitions of that profile within a participant.

Statistical Analysis
The statistical analysis was performed using linear models (LMs) and linear mixed-effect models (LMMs) in R3.6.2 (R Foundation for Statistical Computing, Vienna, Austria). An LM was used to check the constancy of the background activity of all muscles and the torque to exclude potential confounding effects on any of the 49 perturbation profiles. All background scores were standardized within-subject using the Z-score. This standardization was required to avoid heteroscedasticity due to different mV levels of background EMG activity introducing unequal variances across-subjects. An LM with perturbation identifier as fixed effect, that is, using 49 nominal levels, was fit for each background activity separately and the potential effect was evaluated using an ANOVA F test.
LMMs, fitted for each ankle plantarflexor muscle (SOL, GM, and GL) separately, were used to evaluate the dependence of the M1 and M2 stretch reflex. All M1 and M2 responses were normalized within subject by dividing them with the subject mean across all responses, thus expressing M1 and M2 in %EMG mean . This normalization avoided convergence issues due to across-subject differences in response magnitude and variance. A consistent model building strategy was employed across all LMMs to minimize bias within the presented results. First, the fixed effects models were built, which always included an acceleration, velocity and duration predictor to test the main hypotheses. For the M2 response models a two-piece linear predictor for duration was used to fit any nonlinear effects, as observed previously (14,15). The two-piece linear predictor adds a discontinuity to allow the predictor to have a different slope at both sides of this breakpoint (22,23). The breakpoint was placed by minimizing the model residual error using a 5-ms resolution. Such a breakpoint was not added to the M1 model, as we hypothesized that no M1 duration-dependence would appear. In addition, to avoid overfitting, interaction effects were added in full sets per order. Thus, initially all first order interactions, then all second order interactions, etc., as long as model improvements were significant at a = 0.05 using an ANOVA F test. Second, maximum random effects structures were added to the LMMs to allow for between-subject variation of all fixed effects (24). Note, no random effect for intercept was added, because the intercept for each subject was exactly equal to 100%EMG mean due to the applied normalization. The addition of a random effect for every fixed effect induced convergence issues in all models. To achieve convergence, the step-by-step recommendations of Brauer and Curtin (25) were used, selectively removing covariances between random effects as well as any random effect parameters equal to zero. As a result, exact random effect models varied per LMM, for example, SOL M1 model included an acceleration, duration, and acceleration by duration random effect, whereas SOL M2 included an acceleration, velocity, and nonlinear duration random effect.
The main hypotheses about acceleration-, velocity-and duration-dependence of the M1 and M2 responses were evaluated by testing the respective main effects. Conditional main effects, that is, those influenced by interaction effects, were tested across a wide range of conditions to provide insight into the stretch reflex dependencies. For M1, conditional main effects were evaluated at all three velocity and both acceleration levels, as well as the shortest (35-ms) and longest (75-ms) duration. For M2, conditional main effects were evaluated at all three velocity and both acceleration levels as well as the shortest (35-ms), middle (55-ms), and longest (75-ms) duration. The 55-ms duration was added for M2, as the breakpoint for all LMMs was located at this point.
The conditional main effects were tested using a Wald t test with Kenward-Roger correction for DOF and a Bonferroni correction, applied to the P value, for multiple comparison per fixed effect. Unconditional main effects and the interaction effects were tested using a type-II ANOVA F test with Kenward-Roger correction for DOF. Random effects were not used for any statistical inferences and were solely included to improve the DOF and standard error estimates of the fixed effects model.

Simulation Model
A simulation model was implemented (MATLAB 2017 b, MathWorks, Natick, MA) to qualitatively support the experimental results, in a similar fashion to a study by Schuurmans et al. (15). In short, the experimental perturbation profiles were applied to a muscle spindle model to obtain the Ia afferent firing rate. Together with a constant tonic a drive, the Ia firing rate was used as input for a motoneuron pool to simulate neural activity. M1 and M2 response measures were extracted from the motoneuron pool output as in the experimental protocol. The muscle spindle model used within Schuurmans et al.'s study (15) by design lacked an initial burst response after perturbation onset (26). Due to the rapid timing of the stretch reflex response, this initial Ia burst response has been considered as an important contributor to the stretch reflex and especially the M1 response, see Supplemental Fig. S1 (https://doi.org/10.6084/m9.figshare. 13393724) (16,17). Therefore, the muscle spindle model (26) was replaced with a multiscale muscle mechanics model in which this burst does emerge (19).
A multiscale muscle mechanics model was used to obtain Ia afferent firing rate based on applied MTU velocity, and a and fusimotor drive inputs (19). The multiscale element refers to the muscle spindle and muscle-tendon mechanics included within the model. The model has been validated qualitatively, not quantitatively, against well-known muscle spindle firing characteristics in isometric conditions and after ramp-and-hold and triangular stretches (19). The validated model implementation and parameterization was adopted without any changes. The a and fusimotor drives were set to 15% and 70%, respectively, to allow for the initial burst to appear within the Ia afferent response (19).
An integrate-and-fire motoneuron pool model, consisting of 300 neurons, was used to obtain neural output based on Ia firing rate, a drive (42.5 sp/s) and the transport delay (40 ms) (15,20,29). To obtain a suitable model response, the normalized Ia firing rate was scaled (arbitrarily) with a gain of 400 and a drive was set to achieve an approximate background activity of 10 sp/s (15). To serve as input to each fiber of the motoneuron pool, both Ia firing rate and a drive were converted into spike trains via a Poisson process. The model implementation and parameterization were taken directly from studies by Schuurmans et al. and Stienen et al. (15,30) with only a single parameter adaptation. A refractory time constant s r of 5 ms instead of 20 ms was used to better reflect the relative timing of M1 and M2 observed experimentally.
Twelve repetitions of the perturbation profiles were simulated at a 2,048-Hz discrete time frequency to match the experimental protocol. Motoneuron pool output was low-pass filtered (2nd-order, 200 Hz, Butterworth) to smooth the results. M1 and M2 magnitudes were computed using the same data analysis methods as in the experiment. The M1 and M2 analysis windows were placed at 42-57 ms and 57.5-76 ms, respectively, to best accommodate the motoneuron pool output.

RESULTS
We investigated the M1 and M2 stretch reflex response to disentangle previously reported acceleration-, velocity-, and duration-dependence. A total of 49 perturbation profiles were used to elicit stretch reflexes, across 2 acceleration levels, 3 velocity levels, and 10 duration levels. To study our hypotheses, LMMs were fit for M1 and M2 response of the SOL, GM, and GL muscles to these perturbations averaged across 12 repetitions per participant. In addition, we studied an extended set of 167 perturbation profiles within a qualitative simulation environment in support of the experimental findings.

Experiment Reflexive Responses
Participants were able to keep background torque constant at À3 Nm (plantarflexion) throughout the experiment as instructed. Across participants, 4.3% of all stretches was rejected from further analysis, as background torque deviated more than ±0.2 Nm. Per participant, rejection rates varied between 0.5% and 12.1%, similar to rates reported in the study by Schuurmans et al. (15), with a minimum of 7 (of 12) reflex responses used to average across repetitions. For all muscles and the torque, variations in background activity did not consistently differ for any of the 49 perturbation pro- Visual inspection of the time series of the ensemble-averaged SOL reflexive response showed clear effects due to acceleration and duration, see Fig. 3. The time series showed that M1 magnitude increased with acceleration and, contrarily, that M2 magnitude decreased with acceleration. Furthermore, M2 increased with duration for short durations up to around 50 ms. Visual inspection did not show an M1 or M2 velocity-dependence, or M1 duration-dependence. LMMs were used to confirm these observations across all participants, muscles and the entire set of perturbation parameters.

Simulated Reflexive Responses
The stretch reflex arc model allowed for a double burst of activity for both Ia firing rate and motoneuron pool output, see Fig. 4. High acceleration (64 ms, 4 rad/s, 240 rad/s 2 ) showed two Ia firing rate peaks, at 18 and 45 ms, also visible within the bag1 ("dynamic") intrafusal yank profile. A lower acceleration (140 rad/s 2 ) only showed a single peak at 26 ms. The motoneuron pool output showed a double peak output for both accelerations with an M1 response around 49-52 ms and an M2 response around 62.5-69 ms, similar to the experimental results.
Visual inspection of the model time series showed effects of acceleration, velocity and duration on the Ia firing rate and both M1 and M2 magnitude, see Fig. 5. With increased acceleration, the Ia firing rate slope steepened and both peak and steady-state firing rate were reached earlier (Fig. 5A). Moreover, M1 increased with acceleration, whereas M2 showed a nonlinear acceleration dependence (Fig. 5D). With increased velocity, the ascending slope of the Ia firing rate continued to rise longer and toward a higher magnitude, because the perturbation had a longer period of maximum acceleration (Fig. 5B). Both M1 and M2 increased with velocity, although M1 plateaued above 2.0 rad/s (Fig. 5E). Stretch duration only affected the final period of the Ia firing rate with magnitude dropping and reaching steady-state at the set stretch duration (Fig. 5C). Both M1 and M2 increased with duration and reached a plateau value above 23 and 41 ms, respectively (Fig. 5F).
The simulations revealed that the relative timing of the applied perturbation, Ia firing rate, and motoneuron output, as well as motoneuron synchronization were instrumental for the observed dependencies. M1 and M2 were simulated with a 40-ms transport delay and quantified using windows between 42-57 ms and 57.5-76 ms. Therefore, M1 and M2 could only be causally influenced by the perturbation and Ia firing rate between 0-17 ms (M1) and 0-36 ms (M2), see M1/ M2 brackets, Fig. 5, A-C. For example, the Ia firing rate burst around 45 ms observed for high acceleration could not influence either M1 or M2, see Fig. 4 and Fig. 5, A and D. Besides, the plateau observed for the M1 velocity-dependence above 2.0 rad/s could not be explained based on timing (Fig. 5E).
The 2.0 rad/s and 4.0 rad/s perturbations had a different Ia firing rate within the 0-to 17-ms window, see M1 bracket, Fig. 5B. Yet, both M1 magnitudes were equal due to synchronization of firing and refractory periods of all available neurons within the motoneuron pool. Afterward, the increased Ia firing rate within the 0-to 36-ms M2 bracket for the 4.0 rad/s perturbation causes an earlier second synchronized burst (M2) of motoneuron activity with increased magnitude (Fig. 5E).

Short-Latency M1 Dependencies
Experimentally, the increase of SOL M1 magnitude with acceleration was consistently present across participants and perturbations profiles in the LMM, see Fig. 6A and Table 1. The effect size of increasing acceleration ranged from [0.53,0.81] %EMG mean /rad/s 2 (P always <0.001). These differences in effects size were due to the interactions of acceleration with both velocity (F 1,465 = 4.34, P = 0.04) and duration (F 1,465 = 4.57, P = 0.03). The acceleration effect size translated to a modeled difference of 25%EMG mean between the 140 and 175 rad/s 2 levels at 2.5 rad/s and 55 ms. The GM and GL showed similar results, see Supplemental Tables S1-S4, and only results different from the SOL will be highlighted here.
Contrarily, no consistent effects of both velocity (P = [0.37,1]) or duration (P = [0.89,1]) on experimental SOL M1 magnitude were present in the LMM, see Fig. 6, B and C and Table 1. The GL M1 response showed a deviation from the SOL results with an unconditional main effect for duration of À0.10%EMG mean /ms (SE = 0.046) (F 1,10.0 = 5.09, P = 0.05). This duration effect size translated to a modeled difference of only 4.2%EMG mean between the 35-and 75-ms levels.
For the simulation model, M1 dependence showed a split between perturbations below or above the plateau values of 2.0 rad/s and 23 ms, see Fig. 6, D-F and Fig. 8, A-C. Above these 2.0 rad/s and 23 ms threshold values, M1 was unaffected by stretch velocity (Fig. 6E and Fig. 8B) and duration ( Fig. 6F and Fig. 8C) and M1 increased with acceleration ( Fig.  6D and Fig. 8A). The model results matched the experimental dependencies, assuming that the velocity and duration plateau threshold translated to the experimental setting, see  velocity for velocities below 2.0 rad/s (Fig. 8B) and M1 increased with duration for durations below 23 ms (Fig. 8C). Below these threshold values, the acceleration-dependence was limited especially for perturbations with low velocity (Fig. 8A).

Medium-Latency M2 Dependencies
Experimentally, the decrease of SOL M2 with acceleration varied depending on the velocity and duration levels, see Fig. 7A and Table 2. The effect size of increasing acceleration ranged from [À0.72,0.32] %EMG mean /rad/s 2 (P ranged from [<0.001,1]). As for M1, the variation was the result of the interactions of acceleration with both velocity (F 1,448 = 10.6, P = 0.001) and duration (>55 ms: F 1,448 = 3.24, P = 0.07). Specifically, for long durations (>55 ms) and high velocities, (3.0 rad/s) SOL M2 decreased with acceleration. Contrarily, for short durations and low velocities, no effects of acceleration were present. The acceleration effect size translated to a maximum modeled difference of À25%EMG mean between the 140 and 175 rad/s 2 levels at 3.0 rad/s and 75 ms.
The effect of increasing velocity on experimental SOL M2 magnitude depended on the acceleration and duration levels, see Fig. 7B and Table 2. The effect size of increasing velocity ranged from [3.1,69] %EMG mean /rad/s (P ranged from [<0.001,1]). The variation was the result of the interactions with both acceleration, and short and long durations ( 55 ms: F 1,448 = 13.0, P < 0.001; >55 ms: F 1,448 = 6.00, P = 0.01). Specifically, for long durations (>55 ms) SOL M2 increased with velocity, whereas no effects of velocity were present for short durations. The interaction with acceleration did not influence these dependencies. The velocity effect size translated to a maximum modeled difference of 35%EMG mean between the 2.0, 2.5, and 3.0 rad/s levels at 140 rad/s 2 and 75 ms.
To investigate the effect of duration on experimental SOL M2 magnitude a two-piece linear predictor was required, as expected based on results of the upper limb (14,15). The effect size of increasing duration varied due to this nonlinearity and the reported interactions with velocity and acceleration. In general, an increase in SOL M2 with duration was present for short durations ( 55 ms), which leveled off for long durations, see Fig. 7C and Table 2. The positive effect for short durations ranged from [3.0,5.2] %EMG mean /ms (P always <0.001), whereas no effect for long duration was present (P = [0. 16,1]). This translated to a modeled difference of 77%EMG mean between the 35-and 55-ms levels at 2.5 rad/s and 175 rad/s 2 .
The experimental M2 duration-dependence was clearly confirmed within the simulation model, see Fig. 7F and Fig.  8F. The simulated M2 response showed a monotonic increase with duration across all acceleration and velocity levels. Simulated M2 was minimal for the shortest durations, given that most of the Ia afferent response fell within the 0-to 17-ms M1 bracket, see Fig. 5C. Like the experimental results, M2 leveled off for longer durations (41-49 ms) with the exact duration interacting with acceleration and velocity.
The simulated M2 velocity-and acceleration-dependence did not clearly match experimental results, especially within the experimental ranges of 2-3 rad/s and 140-175 rad/s 2 , see Fig. 7, D and E. In general, simulated M2 did show a monotonic increase with velocity across acceleration and duration levels (Fig. 8E). However, in the experimental range the dependence on velocity was limited and sometimes even decreased with velocity. Simulated M2 showed both increasing and decreasing effects with acceleration across the velocity and duration levels ( Fig. 7D and Fig. 8D). However, in the experimental range, simulated M2 mainly increased with acceleration, whereas experimental results showed a steady or decreasing M2. Such a decrease was seen for simulated M2 at higher velocity and accelerations (4.0 rad/s and 175-300 rad/s 2 ). For both velocity and acceleration, simulated dependencies were mainly observed for medium to long durations, as M2 was minimal at short durations.

DISCUSSION
The goal of this paper was to disentangle the duration-, velocity-, and acceleration-dependence and their interactions of the M1 and M2 stretch reflex in the ankle plantarflexors. Experimentally, M1 magnitude increased with acceleration, whereas no effect of velocity or duration was present. These experimental findings were qualitatively replicated with a simulation model for moderate to high velocities and durations. For low velocities or short durations, not included in the experimental protocol, the simulated M1 response did  show velocity-and duration-dependence with a limited acceleration-dependence. Regarding M2, a nonlinear effect of duration was present experimentally as M2 magnitude increased with duration until 55 ms, above which the effect leveled off. M2 magnitude decreased with acceleration and increased with velocity at long durations (>55 ms), whereas no effect of acceleration or velocity was present at short durations ( 55 ms). A monotonic increase in M2 response with duration was replicated with a simulation model. Moreover, the simulation model also showed M2 dependence on acceleration and velocity, although the effect of these dependencies and interaction effects did not clearly match between experiment and simulation.

Short-Latency M1 Dependencies
The M1 response was measured between 49-59 ms experimentally (SOL, median across-participants) and between 42-57 ms in simulation, in line with previously reported latencies (17,31). As such, the stretch perturbation and resulting Ia afferent response could only causally influence the M1 response until 19 ms (experiment) or 17 ms (simulation) after perturbation onset, assuming a 40-ms transport delay (29). The simulation model showed that within this window acceleration, velocity and the shortest durations ( 23 ms) influenced the Ia afferent and M1 responses. Given the previous qualitative validation of the simulation model elements, we expect these results to translate to the experimental setting, although exact results and timings would require detailed parameter optimization (15,19,20). Given the minimum 35-ms duration used experimentally, no M1 durationdependence was expected to be measured. The SOL and GM muscles did indeed not show an experimental duration-dependence, whereas the GL unexpectedly did show a small effect. Overall, causality and timing support the lack of M1 duration-dependence generally observed in experiments, also in previous studies (9,14,15,18).
The observed acceleration-dependence of the M1 response, both experimentally and in simulation, is in line with previous results in the ankle (16,17). The simulation model showed that, through complex biomechanics, increased acceleration caused a steeper slope and higher magnitude of the initial burst in the Ia afferent response. Thus, the M1 acceleration-dependence is linked directly to this initial burst response and resulting motoneuron synchronization as previously simulated and hypothesized (16,17,19).
The observed lack of an M1 velocity-dependence for medium to high velocities, both experimentally and in simulation, contradicts previously published results (8)(9)(10)(11)(12)(13)(14)(15). Importantly, M1 velocity-dependence in the ankle plantarflexors was previously investigated at a larger, mainly higher, range of velocities: 1.5-7.5 rad/s (12) and 1.5-5.0 rad/s (13). Notably, all studies that reported this velocity-dependence did not explicitly keep stretch acceleration constant. Specifically, a scaling of velocity was likely achieved by scaling acceleration, as shown in a study by Dietz et al. (11). This supposed covariation combined with the observed acceleration-dependence in our results suggests that the observed velocity-dependence at medium to high velocities in literature could be explained as actual acceleration effects.
The nonlinear M1 velocity-dependence observed in simulation can be explained through the Ia afferent initial burst response and motoneuron synchronization. The simulation model suggested that, due to timing of all events, M1 velocity-dependence at low velocity was linked to the Ia afferent initial burst with the ascending slope continuing to rise longer toward a higher magnitude. This extends previous results that mainly linked stretch velocity to the Ia afferent dynamic index (19) and shows the importance of simulation models in which this initial burst emerges (15,19,26). Moreover, our simulation results suggested that the M1 velocity-dependence plateau at medium to high velocities was caused through synchronization of motoneuron firing and refractory period (15). As such, the additional excitation due to a higher velocity, at constant acceleration, did not increase M1 magnitude as it fell within a synchronized refractory period of the simulated motoneurons.

Medium-Latency M2 Dependencies
The M2 response was measured between 64-84 ms (SOL, median across participants) experimentally and between 57.5-76 ms in simulation in line with previously reported latencies (17,31). As such, the stretch perturbation and resulting Ia afferent response could only causally influence the M2 response until 44 ms (experiment) or 36 ms (simulation) after perturbation onset, assuming a 40-ms transport delay (29). Both experiment and simulation model showed that within this window acceleration, velocity and duration shape the Ia afferent and M2 responses. The M2 duration-dependence leveled-off at 55 ms (experiment) and 49 ms (simulation). The observed M2 duration-dependence is in line with previous results in the upper limb (9,14,15) and can be explained based on timing and causality, similar to M1. In other words, the M2 response will follow M1 when the stretch duration is applied long enough for a second synchronous burst of activity (M2) to be elicited after the first synchronous burst (M1) (15). The M2 acceleration-and velocity-dependence observed in experiment and simulation roughly matched, albeit at different quantitative values. For both experiment and simulation, the acceleration and velocity dependencies appeared at medium to long durations, as the M2 response was minimal at short durations. The experimental M2 response showed a decrease in magnitude with acceleration at high velocities (3.0 rad/s). This outcome is in line with Finley et al. (17), which showed a nonlinear M2 acceleration-dependence with a decrease between 140 and 175 rad/s 2 . The simulation model qualitatively also showed this nonlinear dependency, however only at a higher velocity (4.0 rad/s) and across a higher acceleration range (175-300 rad/s 2 ).
The experimental M2 response increased with velocity and the simulated M2 response in general also showed a monotonic increase with velocity. In previous studies, contradicting observations were reported either observing an M2 velocity-dependence (8,10,11,13,14) or not (10,12). Moreover, as discussed for M1, stretch acceleration was not explicitly kept constant when changing stretch velocity in previous studies, introducing the acceleration as potential confounder. Overall, the simulation model suggested that the M2 acceleration-and velocity-dependence arise through interaction between Ia afferent initial burst response and the motoneuron dynamics. Specifically, the Ia afferent peak response was often observed around the 17-ms M1 cut-off threshold, thus falling within the synchronized refractory period in the motoneuron pool given the 40-ms transport delay. This mechanism may, for example, explain the M2 decreasing with acceleration as the Ia afferent peak falls earlier resulting in less motoneuron input after the refractory period. The quantitative differences between experiment and simulation are also likely to emerge within this complex physiological interaction. As such, a more detailed investigation of M2 acceleration-and velocity-dependence would be valuable after extensive model parameter optimization.

Stretch Reflex Physiology
The stretch reflex arc consists of causally linked elements, in order: applied stretch perturbation, stretch proprioception, neural transport, and muscle contraction. The simulation model implemented in our study simplifies the arc by using only the Ia afferent as stretch proprioceptor and a monosynaptic motoneuron pool for neural transport. This basic simulation model was able to qualitatively explain most M1 and M2 perturbation dependencies observed experimentally. Essential physiological elements of the simulation model required to achieve these explanations were: MTU and muscle cross-bridge dynamics, Ia afferent sensitivity to intrafusal force and yank, and motoneuron synchronization.
The stretch perturbation at joint scale was translated to the muscle spindle scale through simulation of MTU and muscle cross-bridge mechanics (19). In addition, the resulting intrafusal force and yank profiles, not length, velocity or acceleration, were considered to drive the Ia afferent Model parameters are expressed as %EMGmean per fixed effect unit (±SE). Conditional main effects were tested using a Wald t test with Kenward-Roger correction for degree-of-freedom (DOF) and a Bonferroni correction, applied to the P value, for multiple comparison per fixed effect. Interactions were tested using a type-II ANOVA F test with Kenward-Roger correction for DOF. response (19,32). Combined, these multiscale mechanics determined the Ia afferent initial burst response, which dictated the resulting M1 and M2 responses based on timing. Changes in stretch acceleration and velocity converted to changes in the Ia afferent initial burst, such as modulation of the slope, peak magnitude, and timing, consequently changing M1 and M2.
Motoneuron synchronization of both motoneuron firing and refractory periods translated the single Ia afferent burst into two distinguishable burst within the neural output, M1 and M2 (15). Due to the synchronization, a long enough stretch duration is required to elicit Ia afferent input to trigger the second burst, that is, M2 (15). In addition, the simulation model also suggested that motoneuron synchronization can explain the lack of M1 velocity-dependence for medium to high velocities, as well as observed nonlinearities in M2 acceleration-dependence (17).
Despite general reproduction of experimental results, additional physiological mechanisms could be added to refine the match between experiment and simulation. The monosynaptic Ia afferent pathway alone cannot explain several M2 characteristics (15), such as M2 exceeding M1 (33) or separate modulation of M1 and M2 (34). First, within simulation multiple bursts of Ia afferent activity emerged through the multiscale muscle spindle mechanics, also matching previous experimental observations (35,36). Previous simulation studies with the muscle spindle model also showed these multiple bursts at high stretch accelerations before a more steady-state Ia firing rate is obtained (19). These multiple burst likely emerge due to the cross-bridge cycling kinetics included in the model. Although current timing and causality links M1 and M2 completely to the initial burst, model parameter optimization to experimental results might change relative timing of events (19). Second, additional proprioceptive pathways can contribute to the stretch reflex response. Specifically, multiple studies have shown that the muscle spindle group II afferents are likely to influence the M2 response (18,37). The multiscale muscle spindle model used within this study might offer an interesting framework to also study this type of afferents (19). Third, the M2 response is also likely to originate from a mix of both spinal and transcortical contributions (5, 33).

Study Limitations and Application to Physiological and Clinical Research
The interpretation and generalizability of both experimental and simulation model aspects of our study should be done with care. The experimental perturbation parameter space was limited to 49 unique profiles to avoid muscle fatigue and participant loss of attention. As a result, aspects of M1 and M2 perturbation dependencies were not caught in the experimental data set, as shown through the extended simulation study and other experimental studies, for example, including an extended set of acceleration levels (17). The simulation model was built as a mechanistic model, qualitatively validated with experimental data sets without any additional tuning of the model parameters (15,19). Our basic simulation model misses several known physiological Model parameters are expressed as %EMG mean per fixed effect unit (± SE). Conditional main effects were tested using a Wald t test with Kenward-Roger correction for degree-of-freedom (DOF) and a Bonferroni correction, applied to the P value, for multiple comparison per fixed effect. Interactions were tested using a type-II ANOVA F test with Kenward-Roger correction for DOF. elements and the quantitative fit could be improved through extensive parameter optimization. The combined limitations show in the observed M2 acceleration-and velocity-dependence, which did not fully match between experiment and simulation.
Still, the M1 and M2 dependence on acceleration, velocity, and duration do raise questions beyond the scope of the experiment. Decoupling of perturbation parameters in physiological and clinical research may lead to new or revised conclusions. Most stretch reflex studies at joint level reported on velocitydependence (8)(9)(10)(11)(12)(13)(14)(15), which was potentially confounded by the underlying acceleration-dependence. Similarly, Finley et al. (17) were unaware of any duration-dependency in the ankle plantarflexors, therefore their study design did not control for duration. Given that the duration parameter varied between 50 and 90 ms, thus including duration shorter than 55 ms, their investigation of the M2 acceleration-dependence might have been confounded. At the muscle spindle level, the Ia afferent initial burst and dynamic index characteristics were linked to acceleration and velocity, respectively (19). Applying the systematically designed perturbations to the same model showed that both acceleration and velocity influence the initial burst in distinctive manners.
The clinical evaluation and definitions of spasticity, that is, an exaggerated stretch reflex, have focused on a velocitydependent resistance to passive muscle stretch (2,38). A  Fig. 6 and Fig. 7 is highlighted (gray-shaded area). To enhance visualization, small offsets along the x-axis were used for the individual data points in all subplots. A and D: M1/M2 acceleration-dependence, split across multiple levels of duration (color lightness) and maximum velocity (lines/marker style). B and E: M1/M2 velocity-dependence, split across multiple levels of duration (color lightness) and maximum acceleration (lines/ marker style). C and F: M1/M2 duration-dependence, split across multiple levels of maximum velocity (color lightness), and maximum acceleration (lines/ marker style).
paradigm shift away from a sole velocity-dependence might improve the current understanding of spasticity and its influence on daily living (1). This paradigm shift is further supported by the Ia afferent sensitivity to force and yank (19,32) and successful force-based modeling of spasticity (39). Toward clinical application, the systematic evaluation shows that standardization of perturbation profiles in motorized assessment prototypes is essential. The M1 acceleration-dependence further confirms the hypothesis of Sloot et al. (3,40) that variations in acceleration can account for differences in motorized and manual assessments. Without standardized perturbations, clinical assessments would become device-and protocol-specific, whereas the advantage of adding motorized assessment to clinical practice should lie within its precision and objectivity. This recommendation for standardized tests and consideration of stretch acceleration and duration generalizes to all spasticity evaluation techniques under development, like velocity-dependent stretch reflex thresholds (41) or velocity-based parallel-cascade system identification (42).

Conclusion
Motorized assessment of the stretch reflex or spasticity using ramp-and-hold perturbations should be performed systematically. Experimental protocols should consider all M1 and M2 duration, velocity, and acceleration dependencies, and the interdependence of these perturbation parameters. Using a systematic evaluation, we showed that M1 magnitude depended on stretch acceleration in experiment and simulation. Perturbation parameters outside the experimental scope also showed a nonlinear velocity-and duration-dependence in simulation, explaining the lack of velocity-and duration-dependence observed experimentally. Moreover, we showed that the M2 magnitude in the ankle plantarflexors depended on stretch duration, velocity and acceleration. The simulation model explained these findings using MTU and muscle cross-bridge dynamics, Ia afferent sensitivity to intrafusal force and yank, and motoneuron synchronization. The recommendation for systematic motorized assessment is important for both scientific and clinical applications investigating the physiological origin or effects of neurological disorders on the stretch reflex.

ENDNOTE
At the request of the authors, readers are herein alerted to the fact that additional materials related to this manuscript may be found at https://doi.org/10.4121/c.5571333. These materials are not a part of this manuscript and have not undergone peer review by the American Physiological Society (APS). APS and the journal editors take no responsibility for these materials, for the website address, or for any links to or from it.