Electromyography-driven model-based estimation of ankle torque and stiffness during dynamic joint rotations in perturbed and unperturbed conditions

The simultaneous modulation of joint torque and stiffness enables humans to perform large repertoires of movements, while versatilely adapting to external mechanical demands. Multi-muscle force control is key for joint torque and stiffness modulation. However, the inability to directly measure muscle force in the intact moving human prevents understanding how muscle force causally links to joint torque and stiffness. Joint stiffness is predominantly estimated via joint perturbation-based experiments in combination with system identification techniques. However, these techniques provide joint-level stiffness estimations with no causal link to the underlying muscle forces. Moreover, the need for joint perturbations limits the generalizability and applicability to study natural movements. Here, we present an electromyography (EMG)-driven musculoskeletal modeling framework that can be calibrated to match reference joint torque and stiffness profiles simultaneously via a multi-term objective function. EMG-driven models calibrated on < 2 s of reference torque and stiffness data could blindly estimate reference profiles across 100 s of data not used for calibration. Model calibrations using an objective function comprising torque and stiffness terms always provided less feasible solutions than an objective function comprising solely a torque term, thereby reducing the space of feasible muscle–tendon parameters. Results also showed the proposed framework’s ability to estimate joint stiffness in unperturbed conditions, while capturing differences against stiffness profiles derived during perturbed conditions. The proposed framework may provide new ways for studying causal relationships between muscle force and joint torque and stiffness during movements in interaction with the environment, with broad implications across biomechanics, rehabilitation and robotics.


Introduction
Human movement results from the interaction between the neuromusculoskeletal system and the environment (Winter, 2009). The coordinated activity of all muscles spanning a joint largely defines net joint torque and joint stiffness (Cop et al., 2021), thereby enabling versatile navigation and adaptation to external mechanical demands (Valero-Cuevas, 2016). The ability to determine the muscle force profiles underlying a given movement is crucial to understand how joint torque and stiffness are modulated to enable large repertoires of movements. However, it is currently not possible to directly measure muscle force in the intact moving human in vivo (Herzog, 2017).
Electromyography (EMG)-driven musculoskeletal models (or EMGdriven models) are valuable computational tools to study movement's 2003; Hu et al., 2011). Human joint stiffness is predominantly studied via system identification methods, which require the mechanical perturbation of biological joints (Kearney and Hunter, 1990;van de Ruit et al., 2021). However, the required joint perturbation inherently alters normative neuromuscular function and limits the repertoire of movements, as well as the range of human populations, that can be studied in the first place (Klomp et al., 2013). Moreover, joint-perturbationbased system identification methods provide estimations of joint-level stiffness with no direct link to the underlying muscle forces.
First, we propose a new EMG-driven model that can be calibrated at the joint torque and stiffness levels simultaneously during dynamic ankle joint rotations. We hypothesize this will improve joint stiffness estimation with respect to torque-only calibrated models. Moreover, we hypothesize torque-and-stiffness calibrated models reduce the MTU parameter solution space with respect to torque-only calibrated models, thereby leading to realistic MTU force solutions, i.e., explaining joint stiffness and torque simultaneously. Second, we assess the ability of the proposed modeling framework to estimate joint stiffness in unperturbed conditions, while capturing differences against joint stiffness profiles derived during conventional perturbed conditions.
This study provides a framework to study how muscle force results in joint torque and stiffness modulation during dynamic movements in perturbed and unperturbed conditions. Removing the need for joint perturbation would enable, for the first time, the ability to study the neural control of joint stiffness in natural, unaffected biological systems, while facilitating translation to clinical settings, where joint-perturbation requirements cannot always be met .

Subjects
Five healthy volunteers (age range: 23-30 years, 1 woman) with no self-reported history of neurological or ankle impairments participated in this study. The Natural Sciences and Engineering Sciences Ethics Committee of the University of Twente approved the experimental procedures (reference number: RP 2018-59) and all subjects provided written informed consent. The experiments complied with the Declaration of Helsinki.

Protocol
A single axis dynamometer (MOOG, Nieuw-Vennep, The Netherlands) was used to perturb and measure subjects' right ankle joint angle and torque at 2048 Hz (van 't Veld et al., 2021). The dynamometer's encoder provided an indirect measure of the ankle angle by measuring the position of the footplate. Similarly, a torque sensor placed between the footplate and the dynamometer's actuator served as indirect measure of ankle torque. The right ankle's axis of rotation in the sagittal plane was visually aligned to the actuator's axis of rotation before the start of the experiment. The chair and footplate were adjusted to allow for a knee flexion angle of 45 • and an initial ankle flexion angle of 0 • , i.e., neutral position, respectively. Subjects were instructed to follow a sinusoidal ankle angle target (amplitude: 0.15 rad, frequency: 0.6 Hz) displayed on a monitor, while small pseudo-random angular perturbations (amplitude: 0.03 rad, switching time: 0.15 s) were applied. The dynamometer was controlled via an admittance controller, mimicking a rotational inertia-spring-damper system (inertia: 1 kg m 2 , damping: 2.5 N m s/rad, stiffness: 60 N m/rad). Consequently, ankle torques of approximately ± 10 N m were required to move the dynamometer's footplate through the specified angle profile. Additionally, four subjects also performed the same tracking task, i.e., sinusoidal ankle angle target (amplitude: 0.15 rad, frequency: 0.6 Hz), without being applied any angular perturbations by the dynamometer, i.e., unperturbed trial. Fig. 1 provides an overview of the experimental setup. Fig. 1. Experimental setup. The subject was seated on an adjustable chair with their right foot attached to the dynamometer. The dynamometer perturbed and measured subject's right ankle joint angle and torque during a joint angle tracking task. A monitor provided real-time visual feedback on the angle target (in red) and the subject's ankle angle (in blue). EMGs of five leg muscles were recorded using disc-shaped Ag/AgCl electrodes in a bipolar configuration. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) EMGs were recorded (2048 Hz) using a 32-channel amplifier (Porti, TMSi, Oldenzaal, The Netherlands) and disc-shaped Ag/AgCl electrodes (inter-electrode distance: 24 mm) in a bipolar configuration. Electrodes were placed following SENIAM guidelines (Hermens et al., 2000) on tibialis anterior, soleus, gastrocnemius medialis, gastrocnemius lateralis, and peroneus longus. Three maximum voluntary contraction (MVC) trials of approximately 5 s per muscle group were performed at the start of the experiment to normalize the EMGs. Each MVC trial targeted all recorded muscles and comprised the following tasks: isometric ankle dorsiflexion, plantar flexion, and eversion from a seated position, and standing calf raise.

Data processing
Data processing was performed using MATLAB R2021a (The Mathworks Inc., Natick, MA, USA). EMG signals were band-pass filtered using a zero-lag second-order Butterworth filter (cutoff frequencies: [20 300] Hz), full-wave rectified, and low-pass filtered (cutoff frequency: 3 Hz) using a zero-lag second-order Butterworth filter. The resulting envelopes were normalized by the MVC and resampled at 1024 Hz. In the remainder of the manuscript these normalized EMG envelopes will be referred to as ''muscle excitations''. Measured joint torques and joint angle were low-pass filtered using a zero-lag fourth order Butterworth filter (cutoff frequency: 80 Hz) and resampled at 128 Hz. Processed joint torques and angles will be referred to as reference joint torque and angle profiles.
Reference joint stiffness profiles were obtained applying the Short Data Segments system identification algorithm (Ludvig and Perreault, 2012;Esteban et al., 2019) on the joint torque and angle data measured by the dynamometer. This method uses input and output data across an ensemble of task repetitions, i.e., realizations, to estimate joint stiffness at each time point (van de Ruit et al., 2021). Specifically, the input data were the joint angular displacements due to the external perturbations, and the output data were joint torques measured in response to these perturbations. Reference joint stiffness profiles were filtered using a moving average window of 20 samples.

EMG-driven musculoskeletal modeling
This work extends the Calibrated EMG-Informed Neuromusculoskeletal Modeling Toolbox (CEINMS) we previously developed Durandau et al., 2017). We introduce a new algorithm to calibrate EMG-driven model parameters both at the joint stiffness and torque levels simultaneously. Moreover, we extend previous Hill-type muscle model formulations  to allow for stiffness estimation accounting for MTUs' pennation angle. For an extensive description of the standard EMG-driven modeling formulation via Hilltype muscle models, the reader is referred to Lloyd and Besier (2003), Sartori et al. (2015). The EMG-driven modeling pipeline (Fig. 2) is outlined below. The ''MTU dynamics'' block estimates MTU force and stiffness employing a Hill-type muscle model driven by MTU activation and length with an elastic tendon that uses a Wijngaarden-Dekker-Brent optimization (Brent, 1973) to find the roots of the equilibrium equation between muscle fiber force and tendon force. The ''Joint torque and stiffness computation'' block projects MTU force and stiffness onto the joint level via the MTU moment arms to obtain estimates of joint torque and stiffness. (B) Model calibration: Four parameters per MTU, namely optimal fiber length, tendon slack length, maximum isometric force, and shape factor, are adjusted to best track input reference joint torque and stiffness profiles using the EMG-driven model described in (A). Optimal fiber length and tendon slack length were limited to vary ± 5% of their initial value, the maximum isometric force was scaled with a strength coefficient ∈ (0.3, 2.5) (the MTUs of gastrocnemius lateralis and gastrocnemius medialis, as well as peroneus longus and peroneus brevis, shared the same strength coefficient), and the shape factor could take values ∈ (−3, 0). A simulated annealing optimization routine is used to adjust MTU parameters until the difference between reference (plots in black) and estimated (plots in blue) joint torque and joint stiffness profiles is minimized. The weights of the contributions of joint torque and stiffness, and , respectively, can be tuned to obtain a closer fit to the joint torque or joint stiffness profile. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Activation dynamics. Contrary to previous work (Lloyd and Besier, 2003;Pizzolato et al., 2015), muscle excitations, defined as the normalized EMG envelopes in this study, are mapped into MTU activations ( ) without an intermediate muscle fiber twitch model: where is the experimental muscle excitation, and ∈ (−3, 0) is the shape factor. Five experimental leg muscle excitations were mapped into seven muscle-tendon unit (MTU) activations, i.e., the EMG of the peroneus longus muscle was used to drive the modeled peroneus longus and peroneus brevis MTUs, and the EMG of the tibialis anterior was used to drive the modeled peroneus tertius (i.e., a dorsiflexor) MTU.
MTU kinematics. Joint angles are mapped into MTU length and moment arms using a set of multi-dimensional B-splines (Sartori et al., 2012a).
MTU dynamics. MTU force, , is computed as a function of MTU length, velocity, activation, and pennation angle: where is the muscle's maximum isometric force, (̃), (̃), and (̃) are generic dimensionless active force-length, force-velocity, and passive force-length relationships, respectively,̃and̃are the muscle fiber's normalized length and velocity, respectively, = 0.1 is a damping factor to avoid model singularities when muscles are inactive (Millard et al., 2013), and is the muscle fiber's pennation angle.
MTU stiffness, , is computed as the series arrangement of the equivalent muscle fiber's stiffness in the tendon's line of action, , and the tendon's stiffness, , (Cop et al., 2021): where, based on the work of Jenkins and Bryant (2020), is computed in a generalized way, i.e., considering pennation angle, as: where and are the force and length, respectively, of the muscle fiber along the direction of the tendon's line of action, and are the force and length, respectively, of the muscle fiber along its axis. is computed as: where ( ) is the non-linear tendon force as a function of tendon length , and ( ) is the tendon force as a function of tendon strain ( = ∕ − 1, where is the tendon slack length).
Joint torque and stiffness computation. MTU forces are projected into the joint level to obtain joint torque : where is the moment arm of the th MTU spanning the joint.
The net joint (rotational) stiffness, , is computed as: where represents the stiffness of the th MTU spanning the joint, and is the joint angle.

Model calibration at the joint torque and stiffness levels simultaneously.
The optimal fiber length, tendon slack length, maximum isometric force, and shape factor of each MTU included in the model are calibrated using a simulated annealing optimization routine (Goffe et al., 1994) that minimizes the following multi-term objective function: where , , and are the number of trials, degrees of freedom (DoFs), and samples, respectively, used to calibrate the model, and are weights to the contributions to the objective function of the estimated joint torque and stiffness, respectively, , , and , , are reference joint torque and stiffness values, respectively, , , and , , are modeled joint torque and stiffness values, respectively, and (   ,   ) are the variances of reference joint torque and stiffness profiles, respectively, and is a newly introduced penalty factor that constrains MTUs to operate within a physiological range, i.e., penalizes normalized muscle lengths (̃< 0.5 or̃> 1.5) and negative tendon strains ( < ): where an arbitrary value of = 1000 was used in this study to guaranteẽand were always within the aforementioned boundaries.
For each subject, the generic gait2392 OpenSim model (Delp et al., 1990) was linearly scaled to match their anthropometry. A previously proposed optimization-based method (Modenese et al., 2016) was used to identify initial values for MTU's optimal fiber length and tendon slack length in such a way that their operating range was preserved, with respect to the generic gait2392 OpenSim model, after the linear scaling. Lastly, optimal fiber length, tendon slack length, maximal isometric force and shape factor were further adjusted using our proposed calibration procedure (Section 2.4) to best fit experimental joint torque and stiffness simultaneously (Fig. 2).
For each calibration, only one cycle of the tracking task (≈1.6 s of data) was used. For each subject, calibrations using 35 combinations of and ( and values going from 0 to 1 with a step length of 0.2) were executed. Due to the stochastic component of the simulated annealing algorithm (Goffe et al., 1994), the calibrations using each combination of and were repeated five times to assure convergence to the objective function's global optimum (8). All calibrations were performed on a 64-core processor (AMD Ryzen Threadripper 3990X) and 128 GB RAM workstation, with computation times ranging between 41 and 80 min per calibration. Two calibrated EMG-driven models per subject were selected: the best fit to the torque and stiffness simultaneously, i.e. ''Torque and stiffness'', and the calibration with = 1 and = 0 that best matched the experimental joint torque, i.e. the traditional ''Torque only'' calibration.
The two selected calibrated EMG-driven models were then used to estimate joint torques and stiffness using 100 s (approximately 60 cycles) of new, unseen EMGs and joint angles that were not employed for calibration.

Validation procedures
The results of the simulations and the input data were segmented into cycles and time-normalized between 0 and 100% of the cycle. Three tests were performed.
First test. Validation of the estimated joint torque and stiffness against reference data for the two calibrated EMG-driven models. The rootmean-squared error normalized by the root-mean sum (nRMSE) and the squared Pearson correlation coefficient ( 2 ) were computed to measure similarities in magnitude and shape, respectively (Sartori et al., 2012b).

Second test.
Assessing to what extent a multi-term objective function (8) can identify EMG-driven models that can simultaneously estimate torque and stiffness, compared to the traditional single-term objective function comprising solely differences at the torque level. We computed, for the best fitting calibration of each of the 35 combinations of and , the fitting error at the joint torque level ( = We assessed how many calibrated EMG-driven models obtained joint torque and stiffness fits with root-mean-squared errors (RMSEs) that did not exceed 20% of the RMSE of the best calibration, i.e., the EMGdriven model with the lowest fitting error. Specifically, we compared how many calibrated EMG-driven models obtained acceptable torque fit ( ≤ 1.44⋅min( )) and how many calibrated EMG-driven models obtained acceptable torque and stiffness fit ( + ≤ 1.44⋅min( + )). Additionally, we checked that those calibrated EMG-driven models that obtained a similar torque error underlay different sets of parameters.
Third test. Estimation of joint stiffness via EMG-driven modeling using data from an unperturbed trial and comparison to the results from perturbed counterpart. Using the ''Torque and stiffness'' calibrated EMG-driven model, we qualitatively compared measured joint torques and joint angles from perturbed and unperturbed data. The corresponding joint stiffness estimations were compared via RMSE and 2 metrics. Additionally, a curve analysis was performed using statistical parametric mapping (SPM) (Pataky et al., 2013). Specifically, a paired t-test ( = 0.05) was performed using the spm1d package for MAT-LAB (https://www.spm1d.org/) to identify regions in which the joint stiffness estimations were different with statistical significance. Fig. 3 shows the averaged joint torque and stiffness profiles per subject, derived from EMG-driven models calibrated via the ''Torque only'' and ''Torque and stiffness'' conditions. The joint torque and joint stiffness nRMSEs across all subjects for the ''Torque only'' model ranged between 0.17 and 0.78 (median: 0.38) and 0.19 and 0.92 (median: 0.64), respectively. The joint torque and joint stiffness nRMSEs across all subjects for the ''Torque and stiffness'' model ranged between 0.23 and 0.94 (median: 0.48) and 0.14 and 0.64 (median: 0.32), respectively. The joint torque and joint stiffness 2 values across all subjects for the ''Torque only'' EMG-driven model ranged between 0.56 and 0.97 (median: 0.88) and 3.5 × 10 −5 and 0.92 (median: 0.11), respectively.  The joint torque and joint stiffness 2 values across all subjects for the ''Torque and stiffness'' EMG-driven model ranged between 0.53 and 0.96 (median: 0.83) and 3.4 × 10 −4 and 0.93 (median: 0.49), respectively. Table 1 summarizes nRMSE and 2 values for each subject. Fig. 4 depicts the distribution of nRMSE and 2 values across all subjects and cycles.

Second test
For every subject, the set of and combinations that yielded a similar torque error (i.e., red circles in Fig. 5) always had greater dimensionality than the set of and combinations that yielded a similar total torque + stiffness error + (i.e., blue circles in Fig. 5). Results also showed that all the and combinations that yielded a similar torque error underlay different parameter values and therefore represented actual different model instances (Fig. 6). Across all subjects and modeled MTUs, the median interquartile ranges of the calibrated values of optimal fiber length, tendon slack length, strength coefficient, and shape factor spanned 47%, 42%, 38%, and 50%, respectively, of the permitted values.
From all and combinations with similar torque errors, only a subset of combinations minimized + (≤ 1.44⋅min( + )), i.e. blue Table 1 Normalized root-mean-squared error (nRMSE) and squared Pearson correlation coefficient ( 2 ) values between estimated and reference joint torques and stiffness for each subject. Results expressed as median (interquartile range). Results of the ''Torque only'' calibrated EMG-driven model in gray font, and results of the ''Torque and stiffness'' calibrated EMG-driven model in black font.  Fig. 5. For subject 1, 21 combinations of and resulted in similar torque fits, but only 11 minimized both torque and joint stiffness simultaneously. For subject 2, 23 combinations of and resulted in similar torque fit, but only 9 minimized both torque and joint stiffness simultaneously. For subject 3, 23 combinations of and resulted in similar torque fit, but only 7 minimized both torque and joint stiffness simultaneously. For subject 4, 20 combinations of and resulted in similar torque fit, but only 7 minimized both torque and joint stiffness simultaneously. For subject 5, 19 combinations of and resulted in similar torque fit, but only 17 minimized both torque and joint stiffness simultaneously. Fig. 7 shows, for each subject individually and for the average across all subjects, average joint torque, angle, and stiffness profiles for both the perturbed and unperturbed data. The joint stiffness RMSEs across all subjects between perturbed and unperturbed data ranged between 2.4 × 10 −3 N m/rad and 3.41 N m/rad (median: 0.77 N m/rad). The joint stiffness 2 values across all subjects between perturbed and unperturbed data ranged between 5.4 × 10 −4 and 0.97 (median: 0.67). All subjects, in addition to the average across subjects, showed statistically different joint stiffness profiles in portions of the task cycle: 15%, 30%, 57%, 15% of the whole cycle for subjects 1 -4, respectively, and 46% of the whole cycle for the average across all subjects.

Discussion
We proposed an EMG-driven model-based stiffness estimation methodology, which we validated against the joint perturbation-based Red circles indicate the solutions with a torque error ≤ 1.44⋅min( ), and blue circles indicate the solutions with a total (i.e. torque + stiffness) error + ≤ 1.44⋅min( + ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 6. Box plots, for each modeled muscle-tendon unit (MTU) per subject, of the four muscle parameters that are calibrated: optimal fiber length, tendon slack length, strength coefficient, and shape factor. Optimal fiber length and tendon slack length are normalized with their initial uncalibrated value. Short Data Segments system identification method using data from the same subject population and movement. Differently from joint perturbation-based methods, our approach provides a direct link between joint stiffness and the underlying EMG and muscle force profiles, something crucial to study causalities between neural, muscular, and articular dynamics.
Previous research highlighted the predominant contribution of shortrange stiffness to the model-based estimation of arm end-point stiffness during isometric contractions only (Hu et al., 2011). However, it was never investigated to what extent EMG-driven Hill-type muscle models could estimate reference joint stiffness profiles during dynamic tasks, during which short-range stiffness components would play a less predominant role (Cop et al., 2021).
The first test showed that EMG-driven models calibrated under the ''Torque only'' condition displayed only minimal improvements at estimating reference joint torque profiles, compared to models calibrated under the ''Torque and stiffness'' condition (median nRMSE reduced by 0.1; median 2 increased by 0.05). However, EMG-driven models calibrated under the ''Torque and stiffness'' condition outperformed the ''Torque only'' models at estimating reference joint stiffness profiles (median nRMSE reduced by 0.32; median 2 increased by 0. 38,. Our EMG-driven modeling formulation could estimate reference joint stiffness and torque profiles in a robust way, i.e., EMG-driven models calibrated using 1.6 s of data estimated 100 s of joint torques and stiffness from unseen EMGs and ankle angle profiles ( Table 1). Model-based joint stiffness estimations were in line with literature. Recent work estimated ankle stiffness during dynamometry experiments where subjects tracked a sinusoidal plantar flexion torque while the dynamometer imposed a sinusoidal ankle rotation . Even though the protocol involved higher torque levels, their lowest torque level (i.e., between 0 and 10 N m), yielded a joint stiffness of 25 N m/rad, which matched the stiffness values we obtained at a plantar flexion torque level of 10 N m.
The second test showed there were always more EMG-driven model instances, each of them characterized by different model parameters , that yielded a similar torque fitting error than EMG-driven model instances that yielded a similar total torque + stiffness fitting error (Section 3.2). Moreover, Fig. 5 shows that calibrations only at the torque level, i.e., when ℎ ≠ 0 and = 0, or only at the stiffness level, i.e., ℎ = 0 and ≠ 0, consistently yielded the largest total torque + stiffness fitting error. This provides evidence of the possibility of reducing the MTU parameter solution space by calibrating and validating an EMG-driven model using a multi-term objective function and may facilitate the identification of realistic muscle force solutions that would explain multiple mechanical variables simultaneously.
In the third test, joint stiffness profiles estimated from unperturbed conditions underlay similar trends than stiffness profiles derived from perturbed conditions, with a median RMSE of 0.77 N m/rad and a median 2 of 0.67 (Section 3.3, Fig. 7). This indicates that our EMG-driven model, once calibrated using reference data, was able to estimate realistic joint stiffness profiles without needing to perturb the joint. However, despite similarity between the perturbed and unperturbed conditions, our EMG-driven model was also able to capture subtle differences, enabling, for the first time, the study of joint stiffness in natural, unperturbed conditions. This represents a viable way for understanding joint stiffness modulation during functional movements, e.g., locomotion, where it is not possible to perturb the joints without affecting the underlying neuromechanical processes involved. Moreover, the ability to decode joint stiffness from EMGs and joint angles, without the need to apply external perturbations, might radically change the way wearable assistive robots are myoelectrically controlled.
Follow up studies should extend our proposed methodology to generalize to functional movements. Future work should integrate shortrange stiffness modules that dynamically engage across static and dynamic movements. Future work should systematically investigate what MTU parameters are most sensitive to stiffness (e.g., slopes of the passive force-length curve and the tendon force-strain curve) and enable direct tuning within our proposed calibration method. Previous work explored ankle joint stiffness estimation techniques by combining ultrasonography and system identification during isometric tasks (Jakubowski et al., 2022). Future work will investigate the integration of ultrasonography within our data-driven modeling framework to refine the estimation of MTU states (e.g., activation, length, velocity) and the calibration of parameters at muscle and tendon scales (e.g., tendon slack length) (Dick et al., 2017).
Currently, a limitation of our EMG-driven model's calibration is that it still requires perturbation-based reference joint stiffness profiles for the initial calibration. Nevertheless, the ability of our EMG-driven modeling framework to estimate joint stiffness without perturbing the joints provides a starting point to relax joint-perturbation constraints post-calibration. Moreover, Fig. 3 and Table 1 provide evidence that EMG-driven models (e.g., subjects 1 and 3) could be potentially calibrated in the ''Torque only'' condition, while matching reference joint stiffness profiles. Future work will investigate how to constrain MTU parameters during a ''Torque only'' calibration to enable simultaneous estimation of reference joint stiffness and torque, thereby facilitating our EMG-driven modeling framework's full clinical translation.
Finally, we presented an ankle EMG-driven model that has been calibrated simultaneously at the joint torque and stiffness levels about the plantar-dorsiflexion DoF. Future work will seek to obtain physiologically and dynamically consistent EMG-driven models for joint torque and stiffness estimation across multiple DoFs and joints simultaneously. For this, different joint perturbation and system identification techniques are to be employed (Van der Kooij et al., 2022), and a more extensive calibration comprising a variety of tasks across different joints and planes of movement should be used (Kian et al., 2021). A comprehensive EMG-driven model across all DoFs will enable the implementation of robust and reliable control strategies for assistive devices.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.