Linear Parameter Varying Identification of Dynamic Joint Stiffness during Time-Varying Voluntary Contractions

Dynamic joint stiffness is a dynamic, nonlinear relationship between the position of a joint and the torque acting about it, which can be used to describe the biomechanics of the joint and associated limb(s). This paper models and quantifies changes in ankle dynamic stiffness and its individual elements, intrinsic and reflex stiffness, in healthy human subjects during isometric, time-varying (TV) contractions of the ankle plantarflexor muscles. A subspace, linear parameter varying, parallel-cascade (LPV-PC) algorithm was used to identify the model from measured input position perturbations and output torque data using voluntary torque as the LPV scheduling variable (SV). Monte-Carlo simulations demonstrated that the algorithm is accurate, precise, and robust to colored measurement noise. The algorithm was then used to examine stiffness changes associated with TV isometric contractions. The SV was estimated from the Soleus EMG using a Hammerstein model of EMG-torque dynamics identified from unperturbed trials. The LPV-PC algorithm identified (i) a non-parametric LPV impulse response function (LPV IRF) for intrinsic stiffness and (ii) a LPV-Hammerstein model for reflex stiffness consisting of a LPV static nonlinearity followed by a time-invariant state-space model of reflex dynamics. The results demonstrated that: (a) intrinsic stiffness, in particular ankle elasticity, increased significantly and monotonically with activation level; (b) the gain of the reflex pathway increased from rest to around 10–20% of subject's MVC and then declined; and (c) the reflex dynamics were second order. These findings suggest that in healthy human ankle, reflex stiffness contributes most at low muscle contraction levels, whereas, intrinsic contributions monotonically increase with activation level.


INTRODUCTION
Ankle joint biomechanics can be described by the relationship between the joint position and the torque acting about it, defined as dynamic joint stiffness. It describes the properties of the human actuator and determines (a) the internal load that the central nervous system (CNS) must control and (b) the joint behavior in response to external loads or perturbations. Consequently, a quantitative knowledge of joint stiffness is essential for understanding the normal control of posture and movement and the nature of motor function disorders such as spasticity, rigidity, hypertonia, hypotonia, and flaccidity (Amato and Ponziani, 1999; Bar-On et al., 2014). Also, a good model of joint stiffness is invaluable for the design and control of ankle prostheses and orthoses (Palazzolo et al., 2007).
Joint stiffness modeling has been extensively investigated in the literature (e.g., Kearney et al., 1997;Mirbagheri et al., 2000;Jalaleddini and Kearney, 2011;Sobhani Tehrani et al., 2014). Two distinct physiological mechanisms contribute to joint stiffness: (i) Limb inertia, viscoelasticity of muscle-tendon complex, and active properties of muscle contraction that together define intrinsic stiffness; and (ii) Stretch reflex feedback that changes muscle activation in response to changes in muscle length leading to reflex stiffness. At the human ankle, this has been efficiently modeled with a Parallel-Cascade (PC) structure having separate pathways for intrinsic and reflex stiffness . This study showed that under quasi-stationary conditions, where the joint is perturbed around an operating point (OP) defined by joint position and activation level, the intrinsic stiffness can be modeled by an impulse response function (IRF) and the nonlinear reflex stiffness can be modeled by a Hammerstein system consisting of a static nonlinearity followed by a linear dynamics.
However, numerous quasi-stationary studies, using system identification techniques, demonstrated that both intrinsic and reflex stiffness parameters change drastically and systematically with ankle position and activation level (Weiss et al., 1986;Sinkjaer et al., 1988;Carter et al., 1990;Mirbagheri et al., 2000;Van der Helm et al., 2002;Bar-On et al., 2014;Jalaleddini et al., 2016). Thus, in many functional tasks, like normal gait, where joint position and neural activation continuously change to control movement and counteract external perturbations, joint stiffness will exhibit time-varying (TV) behavior. Furthermore, there is evidence that this TV behavior cannot be predicted simply by interpolating local TI models identified under quasistationary conditions (Kirsch and Kearney, 1997). Therefore, more advanced methodologies are required to identify and characterize joint stiffness during movement or functional tasks.
To this end, a number of approaches have been proposed and used over the years. These include intramuscular mechanism modeling using optimization that minimizes a predefined cost function (Sartori et al., 2015), system identification techniques, or a combination of both (de Vlugt et al., 2010). Methods for identification of TV systems can be divided into four main categories: (i) short segment, (ii) ensemble-based, (iii) timevarying, and (iv) linear parameter varying (LPV).
Short segment methods (Ludvig and Perreault, 2012;Rouse et al., 2014;Jalaleddini et al., 2017) divide non-stationary data into a number of segments with quasi-stationary behavior and identify a time-invariant model for each segment. The segmentation is not always trivial and often requires the TV behavior to be very slow. Ensemble-based methods (MacNeil et al., 1992;Kirsch et al., 1993;Ludvig et al., 2011;Lee and Hogan, 2015) are effective but require many trials with identical TV behavior, which is hard to achieve in many experimental conditions. Moreover, repeating the same task many times may result in fatigue and affect the reliability of estimates. Timevarying identification techniques (Sanyal et al., 2005;Westwick, 2006, 2007;Guarin and Kearney, 2015) use temporal expansion to estimate how the system parameters change continuously with time using data from a single trial; thus simplifying data requirements significantly. However, selecting proper basis functions for temporal expansion is often difficult and the number of model parameters increases significantly if the time-dependent changes are fast; thus reducing the quality of the estimates. Moreover, none of the models identified by these methods can predict the system response to novel trajectories. LPV models have a structure resembling that of linear systems whose parameters change as functions of one or more timedependent signal called scheduling variables (SV). As such, the LPV structure is an excellent candidate for modeling joint stiffness during functional tasks where the TV behavior is mostly due to dependency on neuromuscular variables that vary with time. Also, by relating TV behavior to SVs rather than time, LPV models model the nonlinear mechanisms that generate the TV behavior and thus have the ability to predict the response to novel trajectories. Finally, control theory is well developed for LPV systems (Mohammadpour and Scherer, 2012), which makes LPV models suitable for prostheses and orthoses control.
Despite the significant advantages of LPV models, methods for LPV identification of nonlinear physiological systems have not been studied much. Examples include the LPV modeling of glucose-insulin dynamics in type I diabetes (Cerone et al., 2012) and of the hemodynamic response to profiled hemodialysis (Javed et al., 2010). Our lab has pioneered the use of LPV methods for the identification of joint stiffness. Specifically, Sobhani Tehrani et al. However, these studies were conducted under passive (i.e., at rest) conditions and quantified position dependent changes in stiffness. The study of joint stiffness changes during large timevarying muscle contractions is challenging since neither the muscle activation level nor the voluntary torque are directly measurable as scheduling variable.
In this work, we used the subspace LPV-PC algorithm (Sobhani Tehrani et al., 2014) to characterize changes in both intrinsic and reflex stiffness during isometric, time-varying contractions of the ankle plantarflexors of healthy human subjects. This algorithm, models the intrinsic pathway as a non-parametric LPV impulse response function (LPV IRF) and reflex stiffness as a LPV-Hammerstein cascade of a LPV static nonlinearity and a time invariant (TIV) linear dynamics. The reflex linear dynamic was assumed TIV, similar to previous works (Sinkjaer et al., 1996(Sinkjaer et al., , 1988Ludvig et al., 2011). The scheduling variable, the joint voluntary torque, was estimated from EMG signals using a time-invariant Hammerstein model of EMG-Torque dynamics, which was previously identified using an error-in-variable subspace algorithm. In addition to the experimental examination of the subspace LPV-PC identification method, we also performed Monte-Carlo simulations to demonstrate its accuracy and precision. Figure 1 shows a block diagram of the subspace LPV-PC model with joint angle as input (θ ), total torque as output (TQ tot ), and voluntary torque as scheduling variable (µ). The total torque is the sum of intrinsic (TQ I ), reflex (TQ R ), and voluntary torques (TQ V ), and the colored measurement noise (n). This can be written as:

Problem Formulation
and the stiffness torque is: where, and N represents the total number of samples. The intrinsic stiffness is represented by a LPV IRF model: where h l are the IRF weights that are functions of SV (µ(k)) represented by a basis expansion on the SV: where h ij is the (i, j)-th coefficient for the i-th lag of IRF, g j represents the j-th basis expansion of the SV and n i is the expansion order. Now, rewrite this equation in matrix form to obtain a data equation for the intrinsic pathway; the unknown intrinsic stiffness parameters are: where H l contains the LPV IRF weights for lag l, The basis expansion of the SV can be represented in vector form: and the lagged position inputs with the vector: Then, the input to the intrinsic pathway is constructed by the Kronecker product of Equations (9, 10): Now, rewriting Equation (5) in vector form, the data equation for the intrinsic pathway is: with the regressor: The reflex stiffness is modeled by a differentiator, a delay, and a Hammerstein system comprising a LPV static nonlinearity followed by a time-invariant linear state-space model. The input to the Hammerstein system is the delayed joint velocity (due to reflex delay) denoted by dvel in the equations. The output of Frontiers in Computational Neuroscience | www.frontiersin.org the static nonlinearity is approximated by an orthonormal basis function expansion of the Hammerstein system input, dvel: where, and g i (dvel(k)) is the i-th basis expansion of reflex input (dvel), g j (µ(k)) is the j-th basis expansion of the SV, and ω ij is the coefficient of their products; n p and n r are the expansion orders of the input (dvel) and the SV, respectively. Thus, using basis expansions of the input, the static nonlinearity is converted to n p parallel linear functions, where the expansion weights are dependent on the SV. The vectors of input and SV basis expansions, for reflex pathway, can be written as: with unknown parameters: Thus, the input to reflex linear dynamics becomes: The linear system is modeled using a discrete-time state-space representation of order m: where X(k) is the state vector, z(k) is the input to reflex linear dynamics, and A, B, C, and D are the state-space matrices and: Substituting Equation (17) in Equation (18) yields: where, Combining the data equations for intrinsic and reflex pathways (Equations 12,20), the total joint stiffness can be represented with a Multi-Input-Single-Output (MISO) state-space model: where,  (Verhaegen and Dewilde, 1992) with input and output signals (U T (k) and TQ s (k)) to estimate the order of the system (Equation 22), m. 3. Construct the extended observability matrix using m and the input and output signals, and use it to estimate the state-space matricesÂ R andĈ R . 4. Form the data equation, and isolate the intrinsic and reflex parameters (β I , β R ) in separate terms: where, I and β I are defined in Equations (7, 13), respectively, and: Use orthogonal projection to decompose the total torque into its intrinsic and reflex components: . Use the subspace Hammerstein method described in Sobhani Tehrani et al. (2013b) to estimate the reflex pathway model using dvel(k) as input andTQ R (k) as output.

Methods
We evaluated the performance of the subspace LPV-PC identification algorithm using a simulation study of the LPV-PC model of human's ankle stiffness dynamics (Figure 1). All parameter and nominal values of simulation model were selected based on experimental results reported in literature (Mirbagheri et al., 2000;Jalaleddini et al., 2016).

Model
The intrinsic stiffness was simulated as the LPV IBK model: The inertia (I) and viscosity (B) were set to 0.015 Nm.s 2 /rad and 1.1 Nm.s/rad. The intrinsic elastic parameter (K) and reflex gain and threshold were simulated to have a non-linear behavior with changes in voluntary torque (SV). The linear dynamics of reflex pathway was assumed TIV. Figure 2 demonstrates the simulated parameters. Elasticity was modeled as a polynomial of order 3 for SV. The reflex gain (represented as NL slope in Figure 2C) and threshold (NL threshold, Figure 2D) of reflex Hammerstein system were modeled as polynomial of order 6 for input and a polynomial of order 4 for the SV.
The linear dynamic element of the reflex pathway was assumed to be a second-order low-pass filter with the dynamics: where G = 1 is the system gain, ω n = 25 rad/s is the natural frequency and ζ = 0.9 rad/s is the damping factor. The reflex delay was assumed to be 40 ms. This system was simulated using MATLAB Simulink at 1 kHz for 120 s.

Input and Noise
The input signal was a pseudo random arbitrary level distributed signal (PRALDS) with random switching time uniformly distributed over [250,350] ms, and maximum amplitude equal to 0.05 rad. This input signal was then filtered with a second order Butterworth low-pass filter with cutoff frequency of 30 Hz to represent the actuator dynamics. Output noise was modeled as a white Gaussian signal filtered with a second order Butterworth low-pass filter with cutoff frequency equal to 15 Hz. The noise amplitude was adjusted to produce an average signal-to-noise ratio (SNR) of 10 dB. SNR was calculated as: Figure 3 shows a 4s segment of the position input and noise free and noisy output data, and Figure 4 shows the simulated input (position), scheduling variable (voluntary torque), and torques.

Analysis
To avoid aliasing, all simulation data were filtered with an eighth-order low-pass filter with cutoff frequency of 45 Hz and decimated to 100 Hz before analysis. The intrinsic pathway was identified using a LPV IRF model as described by Equation (5). We calculated the equivalent elasticity of the identified model as the low-frequency (or DC) gain of the LPV IRFs at each SV snapshot. This gain is the steady state value of the integral of identified intrinsic LPV IRF at each SV snapshot.
We assessed the quality of fit by calculating the variance accounted for (VAF): where TQ i represents the noise free simulated torque at time interval i andTQ represented the estimated value; N is the number of samples.
We quantified the quality of identification estimates by using 200 Monte-Carlo trials, each having a new realization of input and noise. The bias and random errors for reflex static nonlinearity estimates were calculated as: where ρ andρ represent true and estimated parameter respectively. Note that both Bias Error and Random Error are also functions of delayed velocity and SV. noise level tested in this simulation study, confirming the efficiency of method in decomposing the total torque into intrinsic and reflex contributions. Figure 7A shows the simulated values of intrinsic pathway elasticity (K) as a function of SV in blue and the mean of 200 Monte-Carlo identification estimates bracketed by two standard deviations of the estimates in red. It is evident that mean of estimates were very close to true value with small variance.

Results
Figures 7B,C show the simulated (blue) and estimated (red) slope and threshold of the estimated nonlinearity extracted from 3D nonlinearity. These values were obtained by finding the best half-wave rectifier (HWR) fit to estimated nonlinearity at each SV using Levenberg-Marquardt method in MATLAB curve fitting toolbox. The red curve shows the mean of 200 Monte-Carlo identification estimates bracketed by two standard deviations of the estimates. The mean of the estimates for slope was very close to simulated values showing that we can accurately retrieve the reflex gain. The estimates of thresholds at some SVs were subject to a maximum of 25% error. There are two explanations for this: (1) the simulated model was different from the identified model, i.e., HWR was simulated and Chebyshev polynomials were used for identification. (2) The distribution of input (velocity for reflex pathway) affects the estimation of threshold. The estimates are expected to be more accurate for an input with rectangular probability distribution. However, these choices were made intentionally in this work to evaluate the performance of the algorithm for a practical case, i.e., true nonlinearity may not be the same as that used for identification for physiological systems, and the actuator dynamics affects the input distribution. Nevertheless, the overall estimated threshold variation trend is very close to the true simulated value. Note that since torque has little power at thresholds, the bias in threshold estimate has little effect on torque prediction.
The LPV nonlinear block of reflex pathway is plotted in Figure 8 in 3D representation; Figure 8A shows the true simulated nonlinearity whereas the average of 200 estimated nonlinear block is plotted in Figure 8B   lower two panels show the bias and random errors for static nonlinearity estimate for 200 simulation trials from top view; both errors were small with maximum bias error occurring around nonlinearity threshold. This is consistent with our estimation of threshold demonstrated in Figure 7C. The maximum bias error was around 10 Nm/rad and the maximum random error was 1 Nm/rad, while the nonlinearity has a maximum gain of 160 Nm/rad. This confirms the efficiency of the proposed algorithm for estimating the LPV static non-linearity.
The frequency response of reflex linear dynamic estimate is demonstrated in Figure 9. The linear system was calculated as a subspace system; the frequency response representation is used for better visualization of accuracy at different frequencies. Both the gain and phase estimates were close to true simulated values.

Methods
The new algorithm was used to characterize the modulation of joint stiffness with activation level in healthy humans performing an isometric torque tracking task of the ankle plantarflexors.
4.1.1. Apparatus Figure 10 shows a schematic of the experimental setup which is described in details in Morier et al. (1990). Subjects lay supine on an experimental table with the left foot attached to a hydraulic actuator using a costume-made fiberglass boot. The neutral position was defined as a 90 degree angle between the foot and shank. Dorsiflexing rotations were taken as positive. The mean ankle angle was set to 0.2 rad.

Subjects
Five healthy subject (one female and four males) aged 26-33 with no history of neuromuscular disorders participated. Subjects gave informed consent to the experimental procedures, which had been reviewed and approved by McGill University Research Ethics Board. Table 1 summarizes the subjects' demographics.

Data Acquisition
EMG signals from tibialis anterior (TA) and triceps surae (TS) including lateral and medial Gastrocnemius muscles were recorded separately using differential surface electrodes. EMGs were amplified and band-pass filtered with a gain of 1,000 and cutoff frequencies 20-2,000 Hz. Ankle torque was low-pass filtered with an eighth-order Bessel filter with cut-off frequency equal to 0.7 Hz in real time and provided to the subject as visual feedback signal. Position, torque and EMG signals were filtered with an anti-aliasing filter at 486.3 Hz, sampled at 1 kHz, and recorded.

Trials
Subjects were instructed to modulate their ankle torque by tracking a visual command signal. The command signal comprised of a sine-wave with a period of 60 s and peak-to-peak amplitude equal to 40% of their maximum voluntary contraction (MVC). Two conditions were examined: 1. Unperturbed trial (UT): a low-amplitude pseudo random binary sequence (PRBS) signal was added to the command signal. No position perturbations were applied. The PRBS perturbation was added to command signal (sine-wave) to provide the rich, persistently excitatory input needed for accurate identification of the EMG-Torque dynamics. 2. Perturbed trial (PT): random perturbations of ankle position were applied by the hydraulic actuator. The perturbation signal was a PRALDS signal with switching rate of 250-350 ms with amplitude of 0.05 rad.
Data were recorded for 120 s at sampling frequency of 1kHz and then decimated to 100 Hz for analysis. Data were examined for evidence of fatigue or co-activation; there was no evidence of either in any of the trial.

Analysis
Identification was performed in three steps: 1. EMG-Torque Dynamics Estimation: We used a time-invariant error-in-variable (EIV) subspace Hammerstein identification algorithm to estimate the dynamic relationship between rectified voluntary Soleus EMG, and torque from UT data. This algorithm provides unbiased estimates of EMG-Torque dynamics in experimental conditions where the feedback is significant as discussed in . This method uses past inputs and outputs as instrumental variables in a manner similar to the subspace Hammerstein identification approach described by Jalaleddini and Kearney (2013 Figure 11 shows the joint position, visual command, full-wave rectified Soleus EMG and measured and predicted voluntary torque from a typical UT trial. The model estimated between Soleus EMG and torque, predicted the torque extremely well; the variance accounted for was 93% for this subject and 92 ± 3% for all subjects. Figure 11E shows the measured and estimated transient torques. These were obtained by filtering the torques with a moving average Butterworth low-pass filter to remove the slow time-varying torques (sine-wave). The VAF for transient response was 81% for this subject. Figure 12 shows the position perturbation, the visual command, and the resulting torque from a typical PT trial. The voluntary torque, estimated from the UT EMG-Torque model is shown in magenta in Figure 12C, superimposed on the total measure torque in blue. The three lower panels show the intrinsic, reflex, and stiffness torques estimated using LPV-PC identification algorithm for the trial segment with largest variation in voluntary torque (SV). Comparing the stiffness torque and that predicted using LPV identification algorithm, it is evident that the LPV method captured the TV behavior of the system well with a VAF of 82% for stiffness torque and 95% for total torque (stiffness + voluntary torque). The total VAF was never <90% in any trial. Figure 13 shows the LPV-PC model estimate for a typical subject. Figure 13A shows the TV behavior of the intrinsic dynamics and how it varies with voluntary torque. Figure 13B shows that the static nonlinearity has a strong uni-directional sensitivity to velocity; the slope varies with voluntary activation increasing from rest to 5 Nm and then decreased. Figures 13C,D show the bode diagram of estimated TIV reflex linear dynamics resembling a second-order low-pass filter system. Figure 14 shows the variation of estimated parameters with voluntary torque for the five subjects. The estimates of intrinsic elasticity and reflex gain (nonlinearity slope) were normalized to their maximum value for the contraction range studied for each subject to allow inter-subject comparison. The original values corresponding to data points in the Figure, can be calculated by multiplying the x-axis value by subject's MVC and y-axis value by their corresponding normalization factor. The MVC and normalization factors for each subject are given in Table 1. The intrinsic elasticity (K) (Figure 14A), monotonically increased with contraction level in all subjects. The reflex gain ( Figure 14B) and threshold (Figure 14C) of the static nonlinearity systematically changed with voluntary contraction. The reflex gain increased with voluntary torque up to 10-30% MVC in different subjects and then decreased. The variation in reflex gain was higher than 50%. The reflex nonlinearity threshold also varied with voluntary torque and was not always zero as assumed in most quasi-stationary studies. Given the results of the simulation study, the estimates of threshold values may be biased but the overall trends are expected to be informative. The reflex linear block was estimated to be a second-order low-pass filter with delay varying between 40 and 45 ms (see Table 1). Figures 14D,E show the gain and phase of reflex linear dynamics represented in frequency domain. The bandwidth of reflex pathway varies between 1.65 and 2.9 Hz in subjects examined in this work.

DISCUSSION
This paper investigated and quantified the effects of voluntary contractions on ankle joint dynamic stiffness and its intrinsic and reflex components. Previous work has demonstrated that voluntary muscle activation causes substantial changes of stiffness during functional tasks (Ludvig and Perreault, 2014). Thus, studying this system during large, continuous variations in voluntary contraction will lead to better understanding of the control of movement. We used a subspace LPV-PC identification algorithm to track stiffness changes during large, isometric voluntary torque contractions. We first validated the method using a Monte-Carlo simulation study. These demonstrated that the method yielded estimates that were accurate, precise (thus reliable) and capable of capturing time-varying stiffness changes similar to those expected from quasi-stationary results, efficiently. We then applied the method to experimental data acquired while healthy human subjects made large, transient voluntary contractions. Our analysis of these data showed that the stiffness dynamics varied significantly with the contraction. We believe that the system identification algorithm used in this study provides an accurate description of intrinsic and reflex stiffness dynamics throughout a voluntary contraction and so can be used to asses the contribution of each pathway to joint mechanics in functional tasks.

Simulation Study
We parameters with torque was obtained by interpolating the results of quasi-stationary experiments with normal human subjects. We used colored output noise with its amplitude adjusted to give an average SNR of 10 dB for each simulation trial. The true experimental noise is expected to be lower than this value (Ludvig et al., 2011). Thus, we evaluated the identification algorithm under condition that is more challenging than that actually seen experimentally. There are two main differences between our simulation study and the experimental conditions: (i) SV estimation: In the simulations we assumed that the voluntary torque could be measured and completely removed from total torque. However, in the experiments, the SV must be estimated from the recoded EMG signal. Any errors in estimating the SV will result in identification performance to be lower than that predicted from the simulations. (ii) Identification model structure: We made two assumptions about the model structure: (1) Stiffness dynamics at the ankle can be represented using a PC model structure; this model has been widely used and shown to be successful in predicting the stiffness torque for both quasi-stationary and TV conditions (Mirbagheri et al., 2000;Sobhani Tehrani et al., 2014;Jalaleddini et al., 2017), (2) The reflex pathway has a delay of 40-45ms; this is shown to be true in many studies (Stein and Kearney, 1995;Kearney et al., 1997;Mirbagheri et al., 2000). There were few assumptions about structures of the elements of the PC model. Thus, for the intrinsic pathway the linear dynamics were modeled as a nonparametric IRF whose length was limited to be less than the reflex delay. For the reflex pathway, the nonlinearity is modeled with an orthonormal expansion whose order minimize the prediction error; the linear dynamics were modeled with a parametric model whose order is determined as part of the identification. The excellent prediction ability of the resulting model demonstrates that it accurately reproduces the observed behavior. It is possible that the true structure is more complex than the PC model (i.e., involve more pathways or have complex pathways such as nonlinear-linear-nonlinear cascade). If so, the model is still useful as an approximation since an arbitrary nonlinear system can be represented by a parallel cascade of block structured elements. However, in such a case, there would no longer be a direct relation between the structure of the model and that of the underlying physiological system; this possibility must be taken into account in the interpretation of the results.

Dynamic Stiffness
Our experimental results showed that stiffness increased with contraction level suggesting that system became more stiff at high contraction levels. The increase in stiffness may be justified by increase in the number of cross-bridges occurring at higher contraction levels. Reflex gain increased going from rest to lowest active level (occurring between 10 and 20% MVC) and then started to decrease. The variation in reflex gain can be explained by recruitment of more muscle fibers at higher contraction levels and existence of an upper-limit in motoneuron pool excitation. The changes in the nonlinearity threshold suggest changes in motoneuron pool excitation threshold with torque levels. These results indicate that contribution of reflex stiffness is highest at low contractions and decreases as contraction level increase, whereas, intrinsic stiffness monotonically increases with contraction level. Note that we did not attempt to parameterize the LPV IRFs for the intrinsic pathway as a second-order system because: (1) intrinsic dynamics may be more complex than the I,B,K model as demonstrated recently in Sobhani et al. (2017) (2) the fitting procedure would involve nonlinear minimization that would introduce an additional source of error. These findings are essential in understanding the role of stretch reflexes during a motor task particularly those involving low contraction levels such as the control of posture and balance. Other works suggested that intrinsic stiffness is not sufficient to maintain stable upright posture (Morasso and Sanguineti, 2002;Moorhouse and Granata, 2007). Our results show that the range of activation where reflex stiffness is significant, varies among subjects and the reflex contribution was substantial in all subjects examined in this study. Comparing our results to those reported in quasi-stationary condition, the reflex maximum contribution was found to occur around 10% MVC and above whereas this was reported to occur at 5% MVC (Mirbagheri et al., 2000). However, it is not clear whether this is due to the dynamics changes due to task or simply because of differences between the subjects who participated in these studies.
In a separate work, we used a similar approach as that used here to estimate the Hammerstein system of reflex pathway, and evaluated the variation in position-reflex EMG dynamics with contraction levels, in isometric condition . It was demonstrated that both gain and threshold of static nonlinearity changed with contraction levels. The results presented in this work combined with that study gives us a comprehensive understanding of how stiffness modulates during isometric TV contractions in plantarflexors of healthy human subjects.
Given the limited dataset required for the subspace LPV-PC identification algorithm used in this study, this can be used toward exploring the effect of some other factors such as contraction history, contraction rate, and contraction trajectory on dynamics of joint stiffness. This can be achieved by repeating the experiment when: (i) the TV torque-matching task starts after a constant activation level is maintained for a short period of time, (ii) use torque-tracking trajectory with different bandwidths (e.g., different periods for sine-wave) as command signal, (iii) use different torque trajectories, e.g., multi-level, and compare the estimated models for each case.

Comparison to Previous Works
The overall trends in our findings agree with the results of quasi-stationary studies. For example, we found that the intrinsic elasticity increased with activation level, similar to the results of Mirbagheri et al. (2000). Also, for reflex gain, We observed a behavior similar to that reported in Jalaleddini et al. (2016). Nonetheless, the magnitudes of the changes were different. We observed 50% increase in stiffness whereas Mirbagheri et al. (2000) reported this to be around 90% for the same range of contraction. Our estimates of reflex gain were similar to those of Mirbagheri et al. (2000), except that we observed a persistence of reflex contribution for a wider range of contraction levels (up to 30% for some subjects). Some other quasi-stationary works reported the maximum reflex contribution to occur around 50% MVC in dorsiflexors (Sinkjaer et al., 1988;Cathers et al., 2004). Based on our experience, this level of activation is very likely to cause fatigue which affects the reliability of results from such experiments. Also, the nominal values reported for maximum reflex contribution based on %MVC might vary among different works due to the differences in measuring the MVCs or the muscle studied. Van Eesbeek et al. (2013) also used the LPV identification to study wrist stiffness in an activation varying task. However, their method was limited to intrinsic estimates and did not decouple the effects of reflex contribution on the total torque variations. Reflex contributions were reported to be minimal in the upper arm (Bennett et al., 1992) but found to be significant in the ankle , wrist (Sinkjaer and Hayashi, 1989), and knee (Ludvig and Perreault, 2014). Consequently, the results of Van Eesbeek et al. (2013) cannot be directly compared to ours. Also, the range of activation is very different in the wrist compared to the ankle. Nevertheless, they showed that the main variation in intrinsic parameters at human wrist was in the elastic parameter, variations in viscosity were small and the inertia was found invariant. This is consistent with our results.
Other studies have used ensemble-based method to evaluate the effect of activation level on joint stiffness. Visser (2010) studied ankle joint stiffness during a sinusoidal torque matching task, where a monotonic increase in elastic parameter with voluntary torque was observed similar to the observation of this study. The main difference with our results was that Visser (2010) found two peaks in the reflex gain at the lowest and highest activation levels. Also, Ludvig and Perreault (2014) used a similar ensemble-based method to study knee stiffness during rapid activation and reported similar results for the elastic parameter. Nonetheless, using ensemble-based methods for activation-varying experiments have a number of drawbacks. It requires the exact same time-varying behavior to be repeated many times while (i) it is extremely difficult to match muscle activation levels between trials, (ii) the muscle recruitment strategy might change to avoid fatigue, (iii) antagonist muscle(s) might be activated in some trials to assist the tracking task, (iv) occurrence of fatigue is inevitable especially if activation levels above 30% are used in the study, (v) the desired torque trajectory needs to be slow enough so that subject can repeat the same task many times, and (vi) system behavior may change from the first experiment to the last one considering the large number of trials required.
The LPV identification algorithm described here, models the underlying dependency of system parameters on torque mean and thus should predict the response to novel trajectories for similar conditions. This predictive ability is a strong asset for studying physiological systems. The experiments described here were not designed to demonstrate this ability but are an important next step. In addition, it is not yet known how this predictive ability depends on the temporal and amplitude properties of the SV. This is an important topic for future work.

Limitations of the study
In this study, we used the subspace LPV-PC algorithm and identified a nonlinear model of both intrinsic and reflex ankle stiffness during isometric, time-varying contractions. The model accurately predicted non-stationary torques recorded from experiments with five healthy subjects. In the identified subspace LPV-PC model, the time-varying behavior of the joint was related to background voluntary torque, instead of time, defined as the scheduling variable. Consequently, it provided insight into functional relationships underlying biomechanics of the joint. Also, the model is expected to predict joint response to novel time trajectories of isometric muscle contractions. However, this study has some limitations too, including: • It assumed that the time-varying behavior of the joint is a function of an a priori known scheduling variable. This assumption was valid for the slow isometric contraction experiments of this study. However, may not hold for other situations such as muscle fatigue, rapid contractions, or neuromuscular disorders where the SV is not well known.
Similarly, it will almost certainly not hold in functional tasks where stiffness parameters depend on multiple variables. For example, during most movements both torque and position change; stiffness parameters are known to depend strongly on both, so it is to be expected that modeling this behavior would require at least two SVs.
• Reflex linear dynamics were assumed to be time-invariant except for its gain that can be modeled by the LPV nonlinearity. This seems to be a valid assumption for healthy subjects performing isometric, slow time-varying contractions (for the contraction range studied in this study) or large imposed movement at rest (Sobhani Tehrani et al., 2014;Jalaleddini et al., 2015). However, it may not be valid for pathological subjects whose reflex dynamics have been shown to change with contraction level (Mirbagheri et al., 2001). Nevertheless, if the subspace LPV-PC identification algorithm is used to analyze a system with TV reflex dynamics, the estimates of intrinsic pathway and corresponding interpretations should remain almost intact. This is because, the subspace LPV-PC identification algorithm uses an orthogonal projection approach to decompose the torque into intrinsic and reflex torques. Thus, any inaccuracy in system structure assumed for reflex dynamics is not expected to affect the estimates of intrinsic dynamics. Rather it would bias estimates of reflex nonlinearity and result in a decrease in torque VAF. Sobhani Tehrani (2017) recently has developed a non-parametric LPV-PC method that can identify SV-dependent changes in reflex dynamics. Future work will use this to investigate the importance of TV changes in reflex dynamics and if this improves the predictions. • The model parameters are assumed to be static functions of the SV while dynamic dependencies may occur in some functional tasks. For the slow isometric contraction trajectory used in this work, the static dependency assumption is expected to be valid. The VAF of its predicted torques supports this assumption. However, assumption must be validated for rapidly changing contractions. In general, if the model parameters depend dynamically on the SV, the LPV identification algorithm would not be expected to predict well. We are not aware of any work investigating potential dynamic dependencies between voluntary torque and joint stiffness parameters. Indeed, the subspace LPV-PC identification algorithm provided the tool needed to investigate such dependencies. • Since the voluntary torque (i.e., the SV) is not directly measurable, we estimated it using an EMG-Torque Hammerstein model, identified from experimental data. The risk is that inaccuracies in the EMG-torque model, and thus the estimated scheduling variable, may bias the identified LPV stiffness model parameters.
Finally, note that this study was performed under open-loop experimental conditions, where the perturbing actuator was many times more stiffness than the ankle. Consequently, the torque generated at the ankle could not change the position of the actuator. This is not the case when subjects interact with compliant loads, where closed-loop conditions may arise. The subspace family of identification algorithms are believed to work with data acquired in closed-loop conditions (Van Wingerden and Verhaegen, 2009); however, validating this with experimental data acquired specifically for LPV-PC modeling of joint stiffness is a subject of future work.

Clinical Significance
The subspace LPV-PC method would be an invaluable tool for objective and quantitative assessment of neuromuscular performance (or impairment) and motor function (or dysfunction). In fact, the early signs of recognizing the clinical benefits of exploiting system identification and modeling approach have recently appeared in the literature (Meskers et al., 2015;Sloot, 2016), where, for example, system identification was used to assess motor dysfunction in children with cerebral palsy. The subspace LPV-method can actually enable and expedite this shift from conventional scoring techniques to model-based clinical assessment, diagnosis, and treatment recommendation. Few of the reasons are: • It works for much more functional tasks compared to quasistationary studies. In addition, the identified LPV model is not just a predictive model. Rather, it provides a coherent representation of the joint biomechanics where the systematic changes are functionally related to variables within the neuromuscular system. • It is far more efficient than the quasi-stationary methods because it requires many fewer trials. For example, in the isometric TV contraction experiment of this study, we used only two trials (UT and PT) to identify the LPV-PC model; whereas the quasi-stationary studies require many more trials to cover the same range of activation levels with a fine resolution. For example, 11 trials are needed to cover activation levels from rest to 40% MVC with a resolution of 2% MVC; thus the LPV method reduces the required number of trials by more than 80%. Such reductions are of utmost importance and value working with patients and in clinical applications. • By estimating the individual elements of the subspace LPV-PC stiffness model, the method distinguishes between the mechanical and reflex contributions to the abnormal joint mechanics, which is very important from a clinical standpoint. Thus, the method will have significant clinical benefits for diagnosis and treatment monitoring of patients suffering from neuromuscular diseases such as cerebral palsy, spinal cord injury, stroke, and Parkinson's disease.

ETHICS STATEMENT
This study was carried out in accordance with the recommendations of McGill University Research Ethics Board with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the McGill University Research Ethics Board.

AUTHOR CONTRIBUTIONS
MAG implemented the simulation study, collected the experimental data, and performed analysis on experimental data. EST developed the identification algorithm. MAG and EST contributed to the execution and drafting of this paper and the work was supervised, reviewed, and approved by REK.

ACKNOWLEDGMENTS
This paper was made possible by NPRP grant #6-463-2-189 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors. This work was also supported by a FQRNT doctorate scholarship to MAG.