==== Front PLoS One PLoS One plos PLOS ONE 1932-6203 Public Library of Science San Francisco, CA USA 10.1371/journal.pone.0282130 PONE-D-23-03468 Research Article Physical Sciences Physics Classical Mechanics Kinematics Physical Sciences Materials Science Material Properties Viscoelasticity Biology and Life Sciences Anatomy Musculoskeletal System Skeleton Skeletal Joints Medicine and Health Sciences Anatomy Musculoskeletal System Skeleton Skeletal Joints Research and Analysis Methods Simulation and Modeling Biology and Life Sciences Biophysics Biophysical Simulations Physical Sciences Physics Biophysics Biophysical Simulations Biology and Life Sciences Computational Biology Biophysical Simulations Biology and Life Sciences Anatomy Body Limbs Legs Medicine and Health Sciences Anatomy Body Limbs Legs Biology and Life Sciences Biomechanics Musculoskeletal Mechanics Biology and Life Sciences Physiology Muscle Physiology Musculoskeletal Mechanics Biology and Life Sciences Anatomy Musculoskeletal System Skeleton Pelvis Hip Medicine and Health Sciences Anatomy Musculoskeletal System Skeleton Pelvis Hip Does joint impedance improve dynamic leg simulations with explicit and implicit solvers? Real-time limb biomechanics https://orcid.org/0000-0002-6141-611X Bahdasariants Serhii Conceptualization Data curation Formal analysis Investigation Methodology Resources Software Validation Visualization Writing – original draft Writing – review & editing 1 Barela Ana Maria Forti Data curation Funding acquisition Investigation Project administration Supervision Writing – review & editing 2 https://orcid.org/0000-0002-6408-9433 Gritsenko Valeriya Data curation Funding acquisition Investigation Software Writing – review & editing 1 3 Bacca Odair Data curation Investigation Resources Writing – review & editing 2 Barela José Angelo Data curation Investigation Writing – review & editing 4 https://orcid.org/0000-0002-5946-6409 Yakovenko Sergiy Conceptualization Data curation Funding acquisition Investigation Methodology Project administration Resources Software Supervision Validation Writing – original draft Writing – review & editing 1 3 5 6 * 1 Department of Human Performance, School of Medicine, West Virginia University, Morgantown, WV, United States of America 2 Institute of Physical Activity and Sport Sciences, Cruzeiro do Sul University, São Paulo, SP, Brazil 3 Department of Neuroscience, School of Medicine, West Virginia University, Morgantown, WV, United States of America 4 Department of Physical Education, São Paulo State University, Rio Claro, SP, Brazil 5 Department of Mechanical and Aerospace Engineering, Benjamin M. Statler College of Engineering and Mineral Resources, West Virginia University, Morgantown, WV, United States of America 6 Department of Chemical and Biomedical Engineering, B.M. Statler College of Engineering and Mineral Resources, West Virginia University, Morgantown, WV, United States of America Naseer Noman Editor Air University, PAKISTAN Competing Interests: The authors have declared that no competing interests exist. * E-mail: seyakovenko@mix.wvu.edu 3 7 2023 2023 18 7 e02821307 2 2023 15 6 2023 © 2023 Bahdasariants et al 2023 Bahdasariants et al https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. The nervous system predicts and executes complex motion of body segments actuated by the coordinated action of muscles. When a stroke or other traumatic injury disrupts neural processing, the impeded behavior has not only kinematic but also kinetic attributes that require interpretation. Biomechanical models could allow medical specialists to observe these dynamic variables and instantaneously diagnose mobility issues that may otherwise remain unnoticed. However, the real-time and subject-specific dynamic computations necessitate the optimization these simulations. In this study, we explored the effects of intrinsic viscoelasticity, choice of numerical integration method, and decrease in sampling frequency on the accuracy and stability of the simulation. The bipedal model with 17 rotational degrees of freedom (DOF)—describing hip, knee, ankle, and standing foot contact—was instrumented with viscoelastic elements with a resting length in the middle of the DOF range of motion. The accumulation of numerical errors was evaluated in dynamic simulations using swing-phase experimental kinematics. The relationship between viscoelasticity, sampling rates, and the integrator type was evaluated. The optimal selection of these three factors resulted in an accurate reconstruction of joint kinematics (err < 1%) and kinetics (err < 5%) with increased simulation time steps. Notably, joint viscoelasticity reduced the integration errors of explicit methods and had minimal to no additional benefit for implicit methods. Gained insights have the potential to improve diagnostic tools and accurize real-time feedback simulations used in the functional recovery of neuromuscular diseases and intuitive control of modern prosthetic solutions. National Institute of Health 1 R03 HD099426-01A1 https://orcid.org/0000-0002-5946-6409 Yakovenko Sergiy FAPESP 2018/04964-8 Barela Ana Maria Forti FAPESP 2018/04964-8 https://orcid.org/0000-0002-5946-6409 Yakovenko Sergiy FAPESP 2018/04964-8 https://orcid.org/0000-0002-6408-9433 Gritsenko Valeriya SY #1 R03 HD099426-01A1, NIH; AB, VG, SY #2018/04964-8, FAPESP. The sponsors played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. Data AvailabilityData and scripts will be published in the public repository on GitHub. Data Availability Data and scripts will be published in the public repository on GitHub. ==== Body pmcIntroduction The central nervous system (CNS) evolved to orchestrate muscular and skeletal actions to produce complex motor behaviors [1]. For instance, in locomotion, intrinsic and sensory feedback signals are integrated with descending visuomotor commands to fine-tune limb stepping and foot placement [2, 3]. To explain this intricate interplay between neural and mechanical systems, previous research has produced a variety of biomechanical models and analyzed their interactions with the environment [4]. For example, a biomechanical model of the arm was shown to explain the directional tuning property of the primary motor cortex neurons, describing why some cortical neurons are most active in particular movement directions [5]. Similar models of neuromechanical interactions are believed to exist within neural pathways in order to overcome transmission delays and nonlinear limb dynamics [6]. The use of embedded dynamical computations for motor control is illustrated in Fig 1. The desired motion transformed through the inverse dynamics yields control signals to muscles (Fig 1A). The oversight over the execution of these motor commands can also be monitored with the use of embedded forward simulations (Fig 1B). Because biomechanical computations can be expressed within neural pathways, models of musculoskeletal transformations are biomimetic and may capture mechanistic cause-effect relationships. 10.1371/journal.pone.0282130.g001 Fig 1 Inverse model of the human motor control (a) and the model validation pipeline (b). From left to right: (1) signal processing—the joint angles θt are found from kinematic equations given the limb positions xt; (2) body dynamics—the joint torques τt are computed using motion equations, given joint angles and their derivatives; (3) muscle dynamics—the muscle forces and computed muscle activations mt are calculated using joint torques (τt,), limb positions (xt), velocities (vt), and musculoskeletal geometry. The model validation process: (1) torques τt, computed using inverse dynamics model (IDM), are compared to reference torques τtref to validate the solution (error: Δτt); (2) joint angles θt*, simulated with forward dynamics model (FDM), are compared to the experimental angles θt to confirm the cause-and-effect relationship between the τt and θt (error: Δθt). Gait abnormalities can be diagnosed and treated more effectively using forward dynamic simulations of body motion [7–9]. However, subject-specific pathomechanics may be required for effective rehabilitation [10, 11]. Simple models may capture the dynamics of the limb in real time, but they may fail to capture the muscle-specific actions that require musculoskeletal morphology [12–18]. Conversely, detailed models with a minimum of 6 degrees of freedom (DOF) and 43 major muscles may require more time to execute [19]. Thus, the engineering challenge in biomechanical modeling is to balance model complexity and computational efficiency. For this reason, we have previously modified Sartori et al. [20] approximation method to derive an objective approximation of muscle dynamics to accurately describe arm and hand muscles in the context of real-time simulations with reduced computational load [21, 22]. This development can be potentially extended to model the muscle dynamics of the leg [23]. However, the fast and accurate simulations of body dynamics remain a challenge [24, 25]. Recently, the use of efficient physical engines [26] has been combined with the musculoskeletal models [19] to develop a new iteration of tools for research, health, and consumer applications [27]. Physical engines use various numerical methods (e.g., Runge-Kutta, implicit/semi-implicit Euler) to integrate ordinary differential equations (ODE) describing multi-DOF body dynamics. Accurate solutions are often achieved by selecting the integration time step (Δt) that satisfies a set tolerance error. Namely, the time step is decreased for rapid system dynamics to maintain integration accuracy and increased for the slow dynamics to speed up the simulation. In real-time decoding, for example, in robot-assisted applications [28], the computing power provides just-in-time or even faster than real-time accurate calculations for movement planning. What is especially challenging is the inclusion of object manipulation in these simulations. Erez et al. compared the performance of real-time engines to demonstrate that each engine is currently optimized for its corresponding use (gaming or robotics) through the difference in system formulation (Cartesian vs. generalized coordinates) and body internal/external forces calculation methods [29]. However, not surprisingly, all simulations fail with the increase in the simulation time step, which indicates a challenge for applications with low sampling rates. Yet, the growing virtual reality use in the post-stroke rehabilitation [30, 31] further necessitates computing body dynamics accurately and swiftly [32]. The use of motion capture systems to acquire body kinematics in these interactive environments is typically performed with fixed sampling frequencies, linked to the low display refresh rate (60 or 120 Hz). Thus, there is a need to improve simulation stability at low sampling rates where numerical inaccuracies result in transient kinematic perturbations. In this study, we developed a biomechanical framework for testing viscoelastic joint constraints and numerical solver types to reduce simulation errors at low sampling rates. We posited that the transient kinematic perturbations due to the large integration time step could be smoothed by adding viscoelastic joint impedance to increase simulation stability. Furthermore, we hypothesized that the viscoelastic compensation depends on the type of a used numerical solver. We tested these hypotheses in a human bipedal model simulating locomotion. We used three numerical solvers with varied combinations of joint impedance to simulate the inverse and forward dynamics of the swing locomotor phase. We determined optimal viscous and elastic pair that stabilized the kinematic solutions and produced close-to-physiological kinetics. In addition, we evaluated feasibility of our approach for real-time applications. Materials & methods This study aimed to test the stabilization effects of viscoelastic compensation, evaluate the performance of numerical solvers, and identify sampling frequencies that produce accurate and stable simulations of the leg dynamics. This was accomplished with the following four steps to develop: (1) subject-specific morphometry scaling of the inertial properties of the body segments, (2) forward and inverse computations for simulating kinematics and kinetics, (3) search for optimal impedance and simulation rate, and (4) validation pipeline for computing numerical errors. Kinematic and kinetic datasets Three datasets were used for the development and validation of the inverse and forward dynamics: (1) experimental marker trajectories, (2) joint angles, and (3) joint torques simulated during the swing locomotor phase. Dataset 1: Marker trajectories The representative three-dimensional marker trajectories were recorded for one healthy individual (male, 26 y/o, 72 kg, 1.75 m). The Institutional Review Board of Cruzeiro do Sul University approved all the procedures, and the participant signed a written consent form prior to the experimental session (institutional review board (CAAE) 02887518.2.0000.8084). A computerized gait analysis system (Vicon, Oxford Metrics, Inc.) with eight infrared cameras tracked reflective markers placed bilaterally on the anterior superior iliac spine, posterior superior iliac spine, medial and lateral epicondyles of the femur, tibialis tuberosity, medial and lateral malleoli, inter-malleolus, and first and fifth metatarsals, following the recommendations of the International Society of Biomechanics [33]. Prior to data acquisition, a T-pose kinematic calibration trial was performed with all reflective markers to register marker clusters relative to the anatomical landmarks. Subsequently, markers were removed from the medial and lateral epicondyles of the femur, tibialis tuberosity, medial and lateral malleoli, and inter-malleolus. The subject was asked to walk at a self-selected comfortable speed without interruption along a 10-m walkway equipped with two force platforms (Kistler, Model 9286BA). Practice was allowed before the data recording session to familiarize the subject with the laboratory environment and the walking condition. From the five recorded trials, 10 unilateral swing phases were extracted using heel-strike and toe-off events, identified based on experimental ground reaction forces and vertical velocity of the foot [34]. Dataset 2: Joint angles The hip, knee, and ankle joint angles, as well as the standing foot rotation with respect to the ground, were calculated from the Dataset 1 using the biomechanics analysis package Visual3D (C-Motion, Inc.). The data were low-pass filtered using a 2nd order Butterworth filter with the cutoff frequency of 6 Hz and resampled with cubic splines [35] at f = 50, 100, 200, 300, 400, 500, 2000 Hz. Dataset 3: Joint torques Joint torques were simulated from the computed joint angles (Dataset 2). Since the numerical solution approaches analytical as the sampling rate increases, torques calculated from angles sampled at 2 kHz were used as a reference (τtref) for all torques simulated at lower frequencies with and without viscoelastic impedance. This impedance was modeled as varied combinations of rotational stiffness and damping coefficients. Model and anatomical scaling A bipedal biomechanical model describing the movement of 17 DOFs was assembled to solve inverse and forward dynamics in Simscape Multibody (R2022b, MathWorks). The proximal to distal segments of the human leg were represented by solid frustums of right circular cones (wide base placed proximally) using subject-specific inertial properties based on validated mathematical relationships [36]. This refinement reduced inertial representation errors to <10% of observed properties. The centers of mass, locations of the hinge points, and inertia were derived for each body segment using Hanavan’s formulations (see S1 Fig). The segment weight was calculated as a fraction of the total human body weight using Winter’s tables [37, p. 86]. The weight of upper body was assigned to the pelvis segment. The axial length of the lower body parts was expressed as a fraction of the subject’s height. The segment geometry was adjusted according to sex and age. Thigh and shank proximal base diameter (d) were calculated from a circumference equation C = πd, where the circumference C was fetched from the sex- and age-based anthropometric reference tables [38, 39]. The proximal base diameter for a foot segment was estimated as a malleolus height. The distal base diameter of all segments was found from the regression equations describing the relationship between axial length and base diameters [36]. The motion of rotational DOFs between the base and follower bodies was expressed as a sequence of time-varying rotational transformations. The knee joint was modeled with a simple 1 DOF revolute primitive. The ankle and hip joints were modeled as a combination of a revolute and universal 2 DOF primitives to avoid the gimbal lock problem (undesired alignment of two perpendicular axes). The foot was rigidly connected to the ground (no slip condition for the leg in stance) close to the metatarsophalangeal (MTP) joint (75% of foot length), and MTP was allowed to rotate in 3 DOFs. This modeled an approximate location of the foot center of pressure during the single stance phase [40]. Since the analysis was focused on the leg in swing, the motion of leg in stance was determined by its kinematics in both forward and inverse simulations. This hybrid design prevented kinetic errors from propagating from the standing to the swinging leg. While 17 DOFs were simulated, we focused our analysis on the seven DOFs of the leg in swing. Body dynamics simulations A built-in recursive algorithm (Simulink R2022b) was used to simulate the inverse dynamics of the modeled body; the explicit 4th-order Runge-Kutta and explicit/implicit Euler methods were used to simulate the forward dynamics. Euler explicit integration is the least accurate and yet the fastest. Euler implicit formulation is better suited to keep the integration error bound to the tolerated error; however, it is slower than the explicit type [41, p. 58]. The accuracy and speed of the Runge-Kutta method are somewhere between those of the two previous methods. For all fixed-step solvers, the integration step is inversely proportional to the sampling rate f. Fig 1B describes the simulation-validation pipeline, where the joint angles (θt, Dataset 2) are calculated from the experimental marker trajectories (Dataset 1) and then used in the inverse dynamics model (IDM) to simulate joint torques (τt, Dataset 3). These torques propagate through the forward dynamics model (FDM) to simulate joint angles (θt*). Simulated and observed angles were compared to assess the accuracy. The kinetic (Δτt) and kinematic (Δθt) errors were calculated as root-mean-square errors (RMSE) between simulated (τt, θt*) and reference / observed signals (τtref, θtref). A dynamic translation technique (Accelerator Mode, Simulink R2022b) allowed us to change input signals without terminating the model execution and increased the simulation efficiency. A standard portable computer running macOS was used for all simulations. Impedance optimization Passive joint impedance was described with the rotational stiffness and damping (k,b) of the spring-dampers used to stabilize multibody simulations [42]. As expected, high impedance not only decreased angular errors but also increased torque errors. This, therefore, presented a sensitivity problem, which involved identifying optimal impedance values to reduce simulation errors. To find the optimal impedance, we searched for the minimal integrated kinematic and kinetic error. Computed kinetic (Δτt) and kinematic (Δθt) errors for each DOF were normalized to peak-to-peak torque and angle ranges, respectively. We covered the grid of multiple k ∈ [0, 10–4, 10–3, 10–2, 10–1] Nm/deg and b ∈ [0, 10–4, 10–3, 10–2, 10–1] Nm/deg/s for multiple experimental swing phases, simulation rates, and numerical solvers: Normalization for each DOF, solver type, and sampling frequency: A=Δθ[N×K×B]/θpp[N×1];T=Δτ[N×K×B]/τpp[N×1], (1) Optimization: (k,b)[D×S×F]=min(A[K×B×D×S×F]+T[K×B×D×S×F]), (2) The superscripts denote dimensionality. N is the number of swing phases; K is the number of stiffness coefficients; B is the number of damping coefficients; D is the number of DOFs; S is the number of solvers; F is the number of simulation rates; θ pp and τ pp are vectors containing joint range of motion and torque peak-to-peak values, respectively; and (k,b) is optimal stiffness and damping matrix. Statistical analysis The normality was tested with the Anderson-Darling test for the following variables: model execution time, kinematic and kinetic errors [43]. The appropriate statistical significance level (α) initially set to 5% was adjusted with the Holm-Bonferroni correction, and the difference between simulations with and without the viscoelastic contribution was tested with the Wilcoxon rank sum test. The linear regression demonstrated the correlations between the execution times and simulation rates. All analyses were performed using the Statistics and Machine Learning Toolbox (MATLAB, R2022b, MathWorks Inc.) Results In this study, we used viscoelastic joint impedance to reduce numerical inaccuracies in simulations of limb dynamics. The impedance, increasing from proximal to distal joints, significantly decreased errors of explicit integration. The beneficial impedance stabilization was negligible when used with the implicit Euler method, which outperformed the explicit Euler and Runge-Kutta methods. The added impedance generated an expected increase in torque errors, however, stabilized kinematics at low sampling rates. Overall, the simulation errors were below 1% of peak-to-peak kinematics, and below 5% of peak-to-peak kinetics. All simulations were executed faster than in real time. Optimal viscoelastic impedance We modified viscoelastic joint impedance to determine the optimal values that minimize both kinetic and kinematic errors in simulations of leg movements. High impedance values increase torque errors (Fig 2A) and decrease kinematic errors (Fig 2B); low impedances have the opposite effect (the results for all DOFs and solver types are demonstrated, S2 and S3 Figs). We, therefore, chose the optimal stiffness and damping sets as generating minimal kinetic and kinematic errors at a given sampling rate. These sets were different for different DOFs and increased in value from proximal to distal joints along the limb, i.e., the viscoelastic contribution was the lowest at the hip and the highest at the ankle. These values are shown with the blue outline in Fig 2 (see S1 and S2 Tables in S1 File contain this information for all DOFs in the model). 10.1371/journal.pone.0282130.g002 Fig 2 Stabilizing impedance affects (a) kinetic and (b) kinematic errors differently across DOFs. The errors are shown for the flexion-extension DOFs of the swinging leg’s hip, knee, and ankle during one representative swing phase simulated with the build-in recursive (a) and implicit Euler (b) methods at 200 Hz. The kinetic performance decreases at high values, and the kinematic performance decreases at low values. The errors obtained with optimal (k,b), specific to DOF, sampling frequency, and solver, are squared in light blue. The improvements in kinematic accuracy with additional impedance were not uniform across three types of numerical solvers. For example, Fig 3A shows no significant difference for kinematic errors in ankle flexion-extension computed with and without impedance across three solvers. However, adding impedance to explicit integration reduced the tendency to numerically "overshoot" exact kinematic solutions in 6 out of 7 analyzed DOFs in the swinging leg. Fig 3B and 3C demonstrate this beneficial improvement in the solution accuracy for the knee flexion-extension and hip abduction-adduction DOFs. The viscoelastic contribution to the implicit Euler method was unnecessary for 5 out of 7 modeled DOFs, with one of these DOFs (internal-external rotation of the hip) benefiting from stabilization only at low sampling frequencies (50–100 Hz). And yet, reinforcing the implicit Euler method with viscoelastic stabilization reduced errors in ankle internal-external rotation and hip adduction-abduction kinematics. Fig 6 shows that the kinetic profiles of these DOFs have frequent sign-switching slopes. These led to higher integration errors and the need for stabilization. Adding impedance to these joints reduced kinematic errors (Fig 3C) but increased kinetic errors. Fig 4A shows that with added impedance, the torque errors increase moderately (marked with *). While the kinetic errors with and without impedance were significantly different in 4 out of 7 examined DOFs, they did not differ at low sampling frequencies in some DOFs. These included ankle flexion-extension (50–100 Hz), its eversion-inversion (50 Hz), and knee flexion extension (50–200 Hz) for which relatively high stiffness (up to 10−1 Nm/deg) or/and damping (up to 10−1 Nm/deg/s) were applied to counter integration inaccuracies in forward simulations (see S1 and S2 Tables in S1 File for details). The viscoelastic contribution was therefore effective for forward simulations with explicit methods and negligible for inverse simulations with fixed low sampling rates. In general, the use of stabilizing impedance was disadvantageous in inverse simulations with high sampling rates. 10.1371/journal.pone.0282130.g003 Fig 3 Stabilizing impedance affects kinematic errors differently across solvers. Demonstrated are three representative cases when impedance had: (a) no stabilizing effect across solvers, (b) no effect for a specific solver, and (c) effect for all solvers. The error distributions are shown for ten swing phases simulated at the sampling frequency of 200 Hz. The bottom and top edges of the box plots show the 25th and 75th percentiles of the error distributions, horizontal lines demonstrate medians, whiskers extend to maximum errors, and outliers are marked with circles. Empty boxes represent distributions computed without impedance, while shaded boxes correspond to errors obtained with optimal impedance (k,b) specific to DOFs and solvers. The horizontal bars indicate significant differences as tested with the Wilcoxon signed-rank test and corrected with the Holm-Bonferroni method (p<1.3×10−3). 10.1371/journal.pone.0282130.g004 Fig 4 Viscoelastic impedance affects inverse (a) and forward (b) dynamic errors across frequencies. The built-in recursive algorithm and three numerical integrators become more accurate with an increasing sampling frequency (empty boxes); viscoelastic impedance increases kinetic and decreases kinematic errors (shaded boxes), as shown for ten swing phases of the knee flexion-extension DOF. For the sampling rates over 200 Hz, kinetic error distributions with and without impedance were significantly different (p<1.83×10−4); the kinematic distributions computed with explicit methods were different across all sampling rates, while the errors obtained with implicit method were the same (p<1.3×10−3). All distributions were compared using the Wilcoxon signed-rank test and corrected with the Holm-Bonferroni procedure. The description of the boxplots is the same as in Fig 3. Abbreviations in legend: EE—explicit Euler, RK—4th order Runge-Kutta, and IE—implicit Euler methods. Numerical solvers The choice of implicit or explicit numerical integration largely determined the accuracy of the forward simulations. The implicit Euler method outperforms explicit methods as demonstrated by low kinematic errors in Figs 3 and 4B. S4 Fig provides further details of the comparative performance for all joints. The best simulation accuracy even at low sampling rates was achieved by the implicit Euler method. Model execution time The forward and inverse dynamic models were executed faster than in real time across all tested sampling frequencies. Model execution times were normalized to the respective swing phase duration. Fig 5 shows the relationship between simulation frequency and execution speed. The real-time speed metric was expressed as the ratio of median values of simulation time (N = 10) to swing durations, i.e., values below 1 are faster than real time. The regression indicates the expected linear increase in execution time with the increase in the sampling rate. The gray dashed lines indicate performance multiples relative to real time. For example, the execution of the inverse model is about 150x and 35x faster than real-time at 50 Hz and 500 Hz, respectively. Due to the difference in algorithmic complexity of integration and differentiation methods, forward simulations are slower than inverse ones across all sampling frequencies. Yet, the forward simulations of the leg model were faster than real time. At 500 Hz, the explicit Euler method (~9x) was faster than the Runge-Kutta (~7x) and implicit Euler (~3.5x) methods. 10.1371/journal.pone.0282130.g005 Fig 5 Model execution speed scales linearly with simulation frequency for (a) inverse and (b) forward dynamic simulations. Built-in recursive solver (shown in red), explicit Euler (EE), Runge-Kutta (RK), and implicit Euler (IE) methods demonstrate faster than real-time performance across tested simulation frequencies. Circles (o) represent the median execution times registered and normalized to the swing phase duration for each simulation. The thin dashed lines indicate performance multiples relative to the real time (when simulation and experimental time are equal). Thick lines illustrate linear regressions; all R-squared values are greater than 0.9887. The equations describing an approximate relationship between the two variables are expressed in legends. Abbreviations in legend are the same as in Fig 4. Kinetic and kinematic trajectories Lower limb dynamics was simulated with forward and inverse models. Fig 6 shows joint torques computed with the built-in recursive solver (A) and respective joint angles calculated with the implicit Euler method at 200 Hz (B). The variability of swing dynamics is depicted with standard deviation from the mean (see shaded areas for ±s.d. in Fig 6). The computed variability is similar to other studies [44–46]. The form and absolute values of the computed average signals are also similar to those demonstrated in previous studies [47–49, 50, p. 165]. We selected a 200 Hz simulation frequency to keep kinematic and kinetic errors within 1% and 5% tolerance, respectively (S4 Fig). Generally, implicit methods do not improve accuracy with the additional impedance. Only explicit methods benefit from it. 10.1371/journal.pone.0282130.g006 Fig 6 Simulated dynamics of the leg in the swing. (a) Kinetic trajectories, simulated at 200 Hz with implicit Euler method and no viscoelastic contribution are shown for ten swing phases. The thick lines correspond to the averages; the shaded areas show ± standard deviation; the thin lines correspond to single trials. (b) Kinematic trajectories; the description is the same as in a. Abbreviations: f-e—flexion-extension; ad-ab—adduction-abduction; int-ext—internal-external rotation; inv-ever—inversion-eversion; in each legend entry, the first movement type corresponds to the negative direction on the vertical axis and the second movement type corresponds to the positive direction. Discussion In this study, we developed a real-time biomechanical framework for testing viscoelastic parameters that reduce the dynamic simulation errors in a human dynamic model of locomotion. The development involved (1) subject-specific body inertia scaling, (2) forward and inverse dynamic computations for simulating kinematics and kinetics, (3) a search for optimal impedance and simulation frequency, and (4) numerical validation and error assessment. Height, weight, age, and sex parameters were used for scaling segment inertia. The simulation errors were small in forward (< 1%) and inverse (< 5%) simulations as computed with implicit Euler and built-in recursive methods, respectively (see also S4 Fig). The inaccuracies in kinematics associated with explicit solvers matched values found in prior research on human locomotion [51, 52]. Added viscoelastic impedance increased the precision of explicit simulations but, in general, decreased the physiological plausibility of the forces computed in inverse simulations. Nevertheless, small errors could be tolerated as they do not surpass the inaccuracies observed in previous studies [7, 53] and result in increased integration accuracy with explicit methods, particularly when the frequency of simulations is low (Fig 5). The implicit integration method outperformed all explicit methods and required no additional compensation from the impedance. Since explicit methods are standard in physics engines, the use of joint impedance could improve their numerical stability; however, implicit methods deliver the best results. This is supported by the comparative demonstration of handling object interactions with standard physics engines and the MuJoCo engine, which utilizes the implicit integration of velocity [26]. Integration method Accurate dynamic simulations of human limbs in real time require an efficient integration method. Analytical solutions provide the best computational efficiency; however, they may not exist for complex mechanical systems requiring fixed- or varied-step numerical methods to simulate body movement [54, p. 5]. Fixed-step methods are typically used to avoid a sample mismatch between input data acquisition and simulation rates; however, they suffer from the accumulating truncation error. Variable-step methods are challenging to implement in real time because they require a predictive extrapolation of real-time inputs and outputs in simulation. The integration methods of approximating ODE can be explicit or implicit. Explicit schemes compute the future state of the system (t+Δt) from its current state (t). However, they are often unstable when simulating ‘stiff’ dynamical systems with slow- and fast-changing dynamical variables [55, p. 614]. For example, the knee swiftly extends during the mid-swing and flexes slowly at the end of the swing of a human gait cycle; the time of ground contact also introduces a sharp transition between swing and stance. Implicit (“backward”) schemes compute approximations using both the current and future states [41, p. 55]. This requires solving for unknowns using a root-finding algorithm (for example, Newton’s method). The additional computational complexity within implicit methods is relatively less taxing than the advantage of accuracy, as shown in Fig 5. We selected the simplest implicit integrator—the implicit Euler method—to solve the forward dynamics of the human leg in swing. Other integrators often used for physical simulations include explicit/semi-implicit Euler and (explicit) 4-th order Runge-Kutta (RK) [56, 57]. The semi-implicit Euler employs a hybrid integration scheme to solve dynamical systems and is an improvement over the least accurate, explicit Euler method. Used in the context of biomechanical simulations, the semi-implicit method can (i) compute joint angles with the fast explicit scheme and (ii) calculate angular velocities with the implicit scheme to stabilize abrupt amplitude changes. Although beneficial for model execution speed, the semi-implicit algorithm may be sensitive to the numerical noise [58]. RK is a typical explicit method and another commonly used integrator. The specific feature of this method is that it solves ODE using a weighted average between four slopes—one at the beginning, two in the middle, and one at the end of the time interval Δt,—improving the overall slope estimate [59, p. 16]. Although decreasing the integration error, RK requires additional computations per iteration, slowing down the simulation. In addition, as all explicit integration schemes, this method performs poorly for stiff ODEs solved with large integration time steps. Our results support the use of the implicit Euler over the explicit Euler and RK methods (Figs 3 and 4). A modern state-of-the-art physical engine MuJoCo uses RK, implicit Euler, and semi-implicit Euler modified to be fully implicit in the velocity integration [26]. While this method was not directly analyzed, our results provide insight into its performance advantages. As we demonstrate in Figs 3 and 6, the solver selection is essential for optimizing a solution to a specific forward dynamic problem. To model human movement, we recommend considering implicit formulations or enforcing the explicit methods with viscoelastic impedance to balance simulation speed, stability, and accuracy (Figs 5, 6 and S2 Fig). Simulation frequency Choosing an optimal simulation frequency is a speed-accuracy tradeoff problem. We solved it by computing inverse-forward dynamics for a range of simulation frequencies and assessing accuracy and execution times. The inverse and forward simulation were performed sequentially, and the numerical errors from both simulations accumulated, as in a similar study [60]. The simulation frequency needs to be minimized with an acceptable level of errors for closed-loop control applications. We chose the error thresholds (<1% kinematic and <5% kinetic) based on our previous evaluation of errors in the approximation of musculoskeletal dynamics [21, 22]. We found forward and inverse computations to be real-time accurate at frequencies as low as 200 Hz (Figs 3 and 5, also see S4 Fig). In simulations of object manipulation, the physics engine MuJoCo could maintain stability even at lower simulation frequencies (62.5 Hz) due to many elegant methodological choices [29]. This suggests further benefits of using hybrid solvers; however, the confounding implementation overhead may be difficult to evaluate in software packages. Viscoelastic impedance Passive mechanical impedance can stabilize physical simulation. Modeled using joint spring-dampers, they generate viscoelastic forces to oppose abrupt changes in position and velocity [61]. However, in inverse simulations, additional impedance may lead to an overestimation of joint forces. We used this tradeoff to constraint our optimization and find the joint impedance that both reduced kinematic errors and minimized force distortion. Our analyses showed that mechanical stabilization did not improve the accuracy of implicitly computed kinematics, but, in general, improved the results of the explicit methods. As a result, the hypothesized increase in simulation stability and the reduction in transient kinematic perturbations were true for the explicit but not for the implicit integrators. It is important to note, however, that the impedance-mediated stabilization is similar to kinematic filtering with the correction to the middle of range of motion set by spring rest length. Here, we applied the standard recommendation for filtering kinematics [37] and observed benefits of this approach. It is likely that in simulations that include abrupt mechanical interactions the impedance-mediated stabilization may augment both implicit and explicit methods. The concept of impedance-mediated stabilization in our study is applied to reducing numerical errors, but it is also theoretically similar to the physiological mechanism. Muscle intrinsic properties and co-contraction generate joint impedance to resist external perturbations. Moreover, the impedance was hypothesized to be the controlled variable for movement and posture control [62–64]. In our previous work, adding expected joint impedance to the musculoskeletal dynamic transformations from joint angles to computed muscle excitations generated activation profiles matching experimentally recorded electromyography [65, p. 63]. While the impedance control has been a controversial theory of motor control in its pure form, the mechanical impedance is an essential feature of muscle coordination. Future studies could extend our implementation to improve simulations of body movement with realistic muscle-driven impedance. Furthermore, this approach may provide an opportunity for improving the state-of-the-art techniques for controlling active myoelectric prosthetics and orthotics used for rehabilitation of movement [66]. Limitations This study was limited to the analysis of the swing phase. The full step analysis including stance is computationally costly typically resolving equations for viscoelastic ground contact with impulse-based velocity-stepping methods, briefly summarized in [29]. Instead, we simulated the standing leg kinematically by attaching it to the ground with a 3-DOF primitive driven by foot segment angles. This simplification introduced kinematic inaccuracies due to the foot roll over the point of contact. The model can be improved in the future by adding translation degrees of freedom and segmenting the foot. Also, the ankle joint was modeled with only rotational DOFs, and the translations could also improve accuracy [67]. Another limitation was the assumption of constant viscoelastic parameters. The time-varying impedance could improve simulation stability beyond the constant values. The magnitude of viscoelastic coefficients could potentially increase to consider muscle co-contraction and further improve stability even at low simulation frequencies. Lastly, we did not simulate forward dynamics with the semi-implicit Euler method because it was not available in Simulink Multibody software. Including this method in the optimal solver search may further improve the speed-accuracy tradeoff for specialized biomechanical models. Conclusions We developed a high-dimensional biomechanical model and an optimization to test viscoelastic stabilization with multiple numerical solver methods and a range of simulation frequencies. We concluded that for the real-time applications, the implicit Euler method with 5 ms time step (200 Hz sampling frequency) and no viscoelastic contribution was optimal. The simulation was kinematically (error <1%) and kinetically (error <5% peak-to-peak) accurate in real time. However, the viscoelastic impedance was essential to achieve accurate forward dynamic solutions with explicit solvers. Supporting information S1 Fig Mathematical model of the human foot inertia. (A) Dimensions and centers of mass of the right circular cones; (B) secondary axes used to derive inertia of the frustum; (C) frustum dimensions used to set up a reference frame; (D) frustum modeling human foot with the reference frame originating at the ankle joint. (PNG) Click here for additional data file. S2 Fig Kinetic (a) and kinematic (b) errors are affected by stabilizing impedance differently. The errors are shown for the hip adduction-abduction, hip internal external rotation, ankle inversion-eversion, and ankle internal external rotations DOFs of the swinging leg during one representative swing phase simulated with the implicit Euler method at 200 Hz. The kinetic performance decreases at high values, and the kinematic performance decreases at low values. The errors obtained with optimal (k,b), specific to DOF, sampling frequency, and solver, are squared in light blue. Abbreviations: ad-ab—adduction-abduction; int-ext—internal-external rotation; inv-ever—inversion-eversion. (PNG) Click here for additional data file. S3 Fig Stabilizing impedance has dissimilar effects on (a,c) kinetic and (b,d) kinematic errors across DOFs and solver types. The errors are shown for the hip flexion-extension, hip adduction-abduction, hip internal-external rotation, knee flexion-extension, ankle flexion-extension, ankle inversion-eversion, and ankle internal-external rotations DOFs of the swinging leg during one representative swing phase simulated with the built-in recursive (a,c), 4th order Runge-Kutta (b), and explicit Euler (d) methods at 200 Hz. The kinetic performance decreases at high values, and the kinematic performance decreases at low values. The errors obtained with optimal (k,b), specific to DOF, sampling frequency, and solver, are squared in light blue. Notice a hundredfold difference in scale between values in Fig 2B and S2B Fig. Abbreviations: f-e—flexion-extension; ad-ab—adduction-abduction; int-ext—internal-external rotation; inv-ever—inversion-eversion. (PNG) Click here for additional data file. S4 Fig Forward (A) and inverse (B) simulations accuracy are shown as a function of sampling rate, viscoelastic contribution, and numerical solver. The simulations in (A) were solved using numerical integrators: EE—explicit Euler method, RK—4th order Runge-Kutta method, and IE—implicit Euler method. Corresponding inverse simulations (B) were solved with a built-in recursive solver—BR. Labels [k,b] marked the error distributions obtained with optimal impedance. (PNG) Click here for additional data file. S1 File (DOCX) Click here for additional data file. We thank Trevor Moon for his insightful critique on model implementation and Matthew Yough for proofreading the article. 10.1371/journal.pone.0282130.r001 Decision Letter 0 Naseer Noman Academic Editor © 2023 Noman Naseer 2023 Noman Naseer https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version0 27 Mar 2023 PONE-D-23-03468 Does joint impedance improve dynamic leg simulations with explicit and implicit solvers?PLOS ONE Dear Dr. Yakovenko, Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process. ============================== ACADEMIC EDITOR: Both reviewers have identifed some minor issues and have made suggestions to improve. The paper can be accepted after minor revisions. ============================== Please submit your revised manuscript by May 11 2023 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file. Please include the following items when submitting your revised manuscript: A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'. A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'. An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'. If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter. If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols. We look forward to receiving your revised manuscript. Kind regards, Noman Naseer, PhD Academic Editor PLOS ONE Journal Requirements: When submitting your revision, we need you to address these additional requirements. 1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at  https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and  https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf 2. Thank you for stating the following financial disclosure:  "SY #1 R03 HD099426-01A1, NIH; AB, VG, SY #2018/04964-8, FAPESP. The sponsors played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript." Please state what role the funders took in the study.  If the funders had no role, please state: "The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript."  If this statement is not correct you must amend it as needed.  Please include this amended Role of Funder statement in your cover letter; we will change the online submission form on your behalf. 3. We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match.  When you resubmit, please ensure that you provide the correct grant numbers for the awards you received for your study in the ‘Funding Information’ section. 4. Thank you for stating the following in the Acknowledgments Section of your manuscript:  "This work was funded by the NIH Award (1 R03 HD099426-01A1) and FAPESP Grant (2018/04964-8). We thank Trevor Moon for his insightful critique on model implementation and Matthew Yough for proofreading the article." We note that you have provided funding information that is not currently declared in your Funding Statement. However, funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.  Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:  "SY #1 R03 HD099426-01A1, NIH; AB, VG, SY #2018/04964-8, FAPESP. The sponsors played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript." Please include your amended statements within your cover letter; we will change the online submission form on your behalf. 5. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide. 6. Please include your full ethics statement in the ‘Methods’ section of your manuscript file. In your statement, please include the full name of the IRB or ethics committee who approved or waived your study, as well as whether or not you obtained informed written or verbal consent. If consent was waived for your study, please include this information in your statement as well.  7. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.  8. Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice. [Note: HTML markup is below. Please do not edit.] Reviewers' comments: Reviewer's Responses to Questions Comments to the Author 1. Is the manuscript technically sound, and do the data support the conclusions? The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented. Reviewer #1: Yes Reviewer #2: Yes **********  2. Has the statistical analysis been performed appropriately and rigorously? Reviewer #1: Yes Reviewer #2: Yes **********  3. Have the authors made all data underlying the findings in their manuscript fully available? The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified. Reviewer #1: Yes Reviewer #2: Yes **********  4. Is the manuscript presented in an intelligible fashion and written in standard English? PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here. Reviewer #1: Yes Reviewer #2: Yes **********  5. Review Comments to the Author Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters) Reviewer #1: The present study is well-written and provides useful information for researchers who are interested in dynamic simulation of human movement. My main concern lies in Section D (lines 185-203). I am a bit confused about the "optimization" problem stated in this section. Did the authors actually solve any optimization problems to find the optimal k and b? Or they just found the minimal integrated kinematic and kinetic error (line 198) from 25 simulations across the ranges of k and b specified in lines 191-192? If the latter is the case, I won't recommend calling it an optimization problem since no actual searching is performed. It seems to me the authors have solved a sensitivity problem not an optimization problem. My specific comments are: 1. Lines 84-85: Shouldn't it be "fail with the increase in the simulation time step"? 2. Eqs (1) and (2): What does each superscript represent in these two equations? For example, what does [NxKxB] mean in Eq. (1)? Please clarify them in the revision. 3. Results section: To help readers better grasp the effects of various choices on the accuracy and stability of the simulation, I would suggest the authors providing a table in the revision to summarize the main findings. 4. Lines 219-220: Have the authors tried to simulate leg swing with the factors that are not optimal? For instance, I will be curious to know the effect of the worst selection of the factors on the accuracy of kinematics and kinetics. Will the resulting kinematics and kinetics errors be much greater than 1% and 5%, respectively? Such information might be useful to researchers having difficulties accessing the optimal choices. 5. Fig.2b: Is angle RMSE actually reported in degree? A RMSE of 5x10^-2 degrees seems very small to me. Just want to confirm it. 6. Fig. 6: Please indicate the sign of each angle in the figure. 7. Fig. S2b: Only two legends are provided here but more than two colors are displayed. Please provide the legends for the other colors. Reviewer #2: The paper proposed a biomechanical framework to test viscoelastic joint constraints and numerical solver types to decrease simulation errors at low sampling rates. Viscoelastic joint impedance to smooth the transient kinematic perturbations due to the large integration time step was used to increase simulation stability. The method was tested on a human bipedal model. The manuscript is well written and structured. However, following suggestion may improve the quality of paper; Discussion section need to be improved. The authors have well described and explained the results but additionally proposed work should be compared with benchmark studies. Moreover English of the manuscript need to be improved as academic writing, and formatting of the paper need to be reviewed as Figure numbering should be corrected, Fig 2A and 2B comes after Fig 3-4 and 5 in text citations. **********  6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files. If you choose “no”, your identity will remain anonymous but your review may still be made public. Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy. Reviewer #1: Yes: Yi-Chung Lin Reviewer #2: Yes: Syed Hammad Nazeer Gilani ********** [NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.] While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step. 10.1371/journal.pone.0282130.r002 Author response to Decision Letter 0 Submission Version1 11 May 2023 Thank you for your comments. We have submitted the point-by-point responses and improved the original submission. Attachment Submitted filename: Response to Reviewers.docx Click here for additional data file. 10.1371/journal.pone.0282130.r003 Decision Letter 1 Naseer Noman Academic Editor © 2023 Noman Naseer 2023 Noman Naseer https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version1 15 Jun 2023 Does joint impedance improve dynamic leg simulations with explicit and implicit solvers? PONE-D-23-03468R1 Dear Dr. Yakovenko, We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements. Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication. An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org. If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org. Kind regards, Noman Naseer, PhD Academic Editor PLOS ONE Additional Editor Comments (optional): Reviewers' comments: Reviewer's Responses to Questions Comments to the Author 1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation. Reviewer #1: All comments have been addressed Reviewer #2: All comments have been addressed ********** 2. Is the manuscript technically sound, and do the data support the conclusions? The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented. Reviewer #1: Yes Reviewer #2: Yes ********** 3. Has the statistical analysis been performed appropriately and rigorously? Reviewer #1: Yes Reviewer #2: Yes ********** 4. Have the authors made all data underlying the findings in their manuscript fully available? The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified. Reviewer #1: Yes Reviewer #2: Yes ********** 5. Is the manuscript presented in an intelligible fashion and written in standard English? PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here. Reviewer #1: Yes Reviewer #2: Yes ********** 6. Review Comments to the Author Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters) Reviewer #1: The authors have adequately addressed all my comments. Nonetheless, there remain a few minor issues that need their attention. For instance, in Figs 5 and 6, instead of using “n.u.”, it should be mentioned in the figure captions. Reviewer #2: The Authors have addressed all of my concerns with the original manuscript. The revised manuscript is ready for publication. ********** 7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files. If you choose “no”, your identity will remain anonymous but your review may still be made public. Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy. Reviewer #1: Yes: Yi-Chung Lin Reviewer #2: Yes: Hammad Nazeer ********** 10.1371/journal.pone.0282130.r004 Acceptance letter Naseer Noman Academic Editor © 2023 Noman Naseer 2023 Noman Naseer https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 23 Jun 2023 PONE-D-23-03468R1 Does joint impedance improve dynamic leg simulations with explicit and implicit solvers? Dear Dr. Yakovenko: I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department. If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org. If we can help with anything else, please email us at plosone@plos.org. Thank you for submitting your work to PLOS ONE and supporting open access. Kind regards, PLOS ONE Editorial Office Staff on behalf of Dr. Noman Naseer Academic Editor PLOS ONE ==== Refs References 1 Wolpert D. M. , “Computational approaches to motor control,” Trends Cogn. Sci., vol. 1 , no. 6 , pp. 209–216, Sep. 1997, doi: 10.1016/S1364-6613(97)01070-X 21223909 2 Wolpert D. M. , Miall R. C. , and Kawato M. , “Internal models in the cerebellum,” Trends Cogn. Sci., vol. 2 , no. 9 , pp. 338–347, Sep. 1998, doi: 10.1016/s1364-6613(98)01221-2 21227230 3 Drew T. , Andujar J.-E. , Lajoie K. , and Yakovenko S. , “Cortical mechanisms involved in visuomotor coordination during precision walking,” Brain Res. Rev., vol. 57 , no. 1 , pp. 199–211, Jan. 2008, doi: 10.1016/j.brainresrev.2007.07.017 17935789 4 Prochazka A. and Yakovenko S. , “The neuromechanical tuning hypothesis,” Prog. Brain Res., vol. 165 , pp. 255–265, 2007.17925251 5 Lillicrap T. P. and Scott S. H. , “Preference Distributions of Primary Motor Cortex Neurons Reflect Control Solutions Optimized for Limb Biomechanics,” Neuron, vol. 77 , no. 1 , pp. 168–179, Jan. 2013, doi: 10.1016/j.neuron.2012.10.041 23312524 6 Franklin D. W. and Wolpert D. M. , “Computational Mechanisms of Sensorimotor Control,” Neuron, vol. 72 , no. 3 , pp. 425–442, Nov. 2011, doi: 10.1016/j.neuron.2011.10.006 22078503 7 De Groote F. et al. , “A physiology based inverse dynamic analysis of human gait: potential and perspectives,” Comput. Methods Biomech. Biomed. Engin., vol. 12 , no. 5 , pp. 563–574, Oct. 2009, doi: 10.1080/10255840902788587 19319704 8 Piazza S. J. , “Muscle-driven forward dynamic simulations for the study of normal and pathological gait,” J. NeuroEngineering Rehabil., vol. 3 , no. 1 , p. 5, Dec. 2006, doi: 10.1186/1743-0003-3-5 16519796 9 Hidler J. and Neckel N. , “Inverse-Dynamics Based Assessment of Gait using a Robotic Orthosis,” in 2006 International Conference of the IEEE Engineering in Medicine and Biology Society, New York, NY: IEEE, Aug. 2006, pp. 185–188. doi: 10.1109/IEMBS.2006.259392 10 Akhundov R. et al. , “Is subject-specific musculoskeletal modelling worth the extra effort or is generic modelling worth the shortcut?,” PLOS ONE, vol. 17 , no. 1 , p. e0262936, Jan. 2022, doi: 10.1371/journal.pone.0262936 35077508 11 Song K. , Anderson A. E. , Weiss J. A. , and Harris M. D. , “Musculoskeletal models with generic and subject-specific geometry estimate different joint biomechanics in dysplastic hips,” Comput. Methods Biomech. Biomed. Engin., vol. 22 , no. 3 , pp. 259–270, Feb. 2019, doi: 10.1080/10255842.2018.1550577 30663342 12 Mochon S. and McMahon T. A. , “Ballistic walking,” J. Biomech., vol. 13 , no. 1 , pp. 49–57, Jan. 1980, doi: 10.1016/0021-9290(80)90007-X 7354094 13 Onyshko S. and Winter D. A. , “A mathematical model for the dynamics of human locomotion,” J. Biomech., vol. 13 , no. 4 , pp. 361–368, Jan. 1980, doi: 10.1016/0021-9290(80)90016-0 7400165 14 Mena D. , Mansour J. M. , and Simon S. R. , “Analysis and synthesis of human swing leg motion during gait and its clinical applications,” J. Biomech., vol. 14 , no. 12 , pp. 823–832, 1981, doi: 10.1016/0021-9290(81)90010-5 7328089 15 Pandy M. G. and Berme N. , “Synthesis of human walking: a planar model for single support,” J. Biomech., vol. 21 , no. 12 , pp. 1053–1060, 1988, doi: 10.1016/0021-9290(88)90251-5 2577951 16 Taga G. , “A model of the neuro-musculo-skeletal system for human locomotion,” Biol. Cybern., vol. 73 , no. 2 , pp. 97–111, 1995.7662771 17 Gilchrist L. A. and Winter D. A. , “A multisegment computer simulation of normal human gait,” IEEE Trans. Rehabil. Eng., vol. 5 , no. 4 , pp. 290–299, Dec. 1997, doi: 10.1109/86.650281 9422454 18 Full R. J. and Koditschek D. E. , “Templates and anchors: neuromechanical hypotheses of legged locomotion on land,” J. Exp. Biol., vol. 202 , no. 23 , pp. 3325–3332, 1999.10562515 19 Delp S. L. et al. , “OpenSim: open-source software to create and analyze dynamic simulations of movement,” IEEE Trans. Biomed. Eng., vol. 54 , no. 11 , pp. 1940–1950, 2007.18018689 20 Sartori M. , Reggiani M. , van den Bogert A. J. , and Lloyd D. G. , “Estimation of musculotendon kinematics in large musculoskeletal models using multidimensional B-splines,” J. Biomech., vol. 45 , no. 3 , pp. 595–601, Feb. 2012, doi: 10.1016/j.jbiomech.2011.10.040 22176708 21 Sobinov A. , Boots M. T. , Gritsenko V. , Fisher L. E. , Gaunt R. A. , and Yakovenko S. , “Approximating complex musculoskeletal biomechanics using multidimensional autogenerating polynomials,” PLOS Comput. Biol., vol. 16 , no. 12 , p. e1008350, Dec. 2020, doi: 10.1371/journal.pcbi.1008350 33326417 22 Smirnov Y. , Smirnov D. , Popov A. , and Yakovenko S. , “Solving musculoskeletal biomechanics with machine learning,” PeerJ Comput. Sci., vol. 7 , p. e663, Aug. 2021, doi: 10.7717/peerj-cs.663 34541309 23 Modenese L. and Kohout J. , “Automated Generation of Three-Dimensional Complex Muscle Geometries for Use in Personalised Musculoskeletal Models,” Ann. Biomed. Eng., vol. 48 , no. 6 , pp. 1793–1804, Jun. 2020, doi: 10.1007/s10439-020-02490-4 32185569 24 Boots M. T. , “A Biomimetic Approach to Controlling Restorative Robotics,” PhD, West Virginia University Libraries, 2019. doi: 10.33915/etd.7384 25 Manukian M. , “Real-time simulation of arm and hand dynamics using ANN,” p. 53, 2021. 26 Todorov E. , Erez T. , and Tassa Y. , “MuJoCo: A physics engine for model-based control,” in 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, Vilamoura-Algarve, Portugal: IEEE, Oct. 2012, pp. 5026–5033. doi: 10.1109/IROS.2012.6386109 27 Caggiano V. , Wang H. , Durandau G. , Sartori M. , and Kumar V. , “MyoSuite—A contact-rich simulation suite for musculoskeletal motor control,” 2022, doi: 10.48550/ARXIV.2205.13600 28 Abdullah H. A. , Tarry C. , Datta R. , Mittal G. S. , and Abderrahim M. , “Dynamic biomechanical model for assessing and monitoring robot-assisted upper-limb therapy,” J. Rehabil. Res. Dev., vol. 44 , no. 1 , pp. 43–61, Jan. 2007, doi: 10.1682/JRRD.2006.03.0025 17551857 29 Erez T. , Tassa Y. , and Todorov E. , “Simulation tools for model-based robotics: Comparison of Bullet, Havok, MuJoCo, ODE and PhysX,” in 2015 IEEE International Conference on Robotics and Automation (ICRA), May 2015, pp. 4397–4404. doi: 10.1109/ICRA.2015.7139807 30 Darekar A. , McFadyen B. J. , Lamontagne A. , and Fung J. , “Efficacy of virtual reality-based intervention on balance and mobility disorders post-stroke: a scoping review.,” J. Neuroengineering Rehabil., vol. 12 , no. 1 , p. 46, 2015, doi: 10.1186/s12984-015-0035-3 25957577 31 Laver K. E. , Lange B. , George S. , Deutsch J. E. , Saposnik G. , and Crotty M. , “Virtual reality for stroke rehabilitation,” Cochrane Database Syst. Rev., vol. 2018 , no. 1 , Nov. 2017, doi: 10.1002/14651858.CD008349.pub4 29156493 32 Torrisi M. et al. , “Beyond motor recovery after stroke: The role of hand robotic rehabilitation plus virtual reality in improving cognitive function,” J. Clin. Neurosci., vol. 92 , pp. 11–16, Oct. 2021, doi: 10.1016/j.jocn.2021.07.053 34509235 33 Wu G. et al. , “ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion—part I: ankle, hip, and spine,” J. Biomech., vol. 35 , no. 4 , pp. 543–548, Apr. 2002, doi: 10.1016/S0021-9290(01)00222-6 11934426 34 O’Connor C. M. , Thorpe S. K. , O’Malley M. J. , and Vaughan C. L. , “Automatic detection of gait events using kinematic data,” Gait Posture, vol. 25 , no. 3 , pp. 469–474, Mar. 2007, doi: 10.1016/j.gaitpost.2006.05.016 16876414 35 Ga W. , “Data smoothing and differentiation procedures in biomechanics,” Exerc. Sport Sci. Rev., vol. 10 , 1982, Accessed: Sep. 15, 2022. [Online]. Available: https://pubmed.ncbi.nlm.nih.gov/6749517/ 36 Hanavan E. P. Jr , “A mathematical model of the human body,” Air Force Aerospace Medical Research Lab Wright-patterson AFB OH, 1964. 37 Winter D. A. , Biomechanics and motor control of human movement. John Wiley & Sons, 2009. 38 Fryar C. D. , Carroll M. D. , Gu Q. , Afful J. , and Ogden C. L. , “Anthropometric reference data for children and adults: United States, 2015–2018,” 2021. 39 McDowell M. A. , Fryar C. D. , Hirsch R. , and Ogden C. L. , “Anthropometric reference data for children and adults: U.S. population, 1999–2002,” Adv. Data, no. 361, pp. 1–5, Jul. 2005. 40 Fuchioka S. , Iwata A. , Higuchi Y. , Miyake M. , Kanda S. , and Nishiyama T. , “The Forward Velocity of the Center of Pressure in the Midfoot is a Major Predictor of Gait Speed in Older Adults,” Int. J. Gerontol., vol. 9 , no. 2 , pp. 119–122, Jun. 2015, doi: 10.1016/j.ijge.2015.05.010 41 Kreiss H. and Ortiz O. E. , Introduction to numerical methods for time dependent differential equations. Hoboken, New Jersey: Wiley, 2014. 42 Andrews S. , Teichmann M. , and Kry P. G. , “Geometric Stiffness for Real-time Constrained Multibody Dynamics,” Comput. Graph. Forum, vol. 36 , no. 2 , pp. 235–246, 2017, doi: 10.1111/cgf.13122 43 Stephens M. A. , “EDF Statistics for Goodness of Fit and Some Comparisons,” J. Am. Stat. Assoc., vol. 69 , no. 347 , pp. 730–737, Sep. 1974, doi: 10.1080/01621459.1974.10480196 44 Roos P. E. and Dingwell J. B. , “Influence of simulated neuromuscular noise on movement variability and fall risk in a 3D dynamic walking model,” J. Biomech., vol. 43 , no. 15 , pp. 2929–2935, Nov. 2010, doi: 10.1016/j.jbiomech.2010.07.008 20708189 45 Winter D. A. , “Kinematic and kinetic patterns in human gait: Variability and compensating effects,” Hum. Mov. Sci., vol. 3 , no. 1 , pp. 51–76, Mar. 1984, doi: 10.1016/0167-9457(84)90005-8 46 Mitchell A. and Martin A. E. , “Quantifying the effect of sagittal plane joint angle variability on bipedal fall risk,” PLOS ONE, vol. 17 , no. 1 , p. e0262749, Jan. 2022, doi: 10.1371/journal.pone.0262749 35081142 47 Grimmer M. , Eslamy M. , and Seyfarth A. , “Energetic and Peak Power Advantages of Series Elastic Actuators in an Actuated Prosthetic Leg for Walking and Running,” Actuators, vol. 3 , pp. 1–19, Feb. 2014, doi: 10.3390/act3010001 48 Millard M. , Kubica E. , and McPhee J. , “Forward dynamic human gait simulation using a SLIP target model,” Procedia IUTAM, vol. 2 , pp. 142–157, Dec. 2011, doi: 10.1016/j.piutam.2011.04.015 49 Font-Llagunes J. , Pàmies-Vilà R. , Alonso J. , and Lugrís U. , “Simulation and design of an active orthosis for an incomplete spinal cord injured subject,” Procedia IUTAM, vol. 2 , pp. 68–81, Jul. 2011, doi: 10.1016/j.piutam.2011.04.007 50 Robertson D. G. E. , Caldwell G. E. , Hamill J. , Kamen G. , and Whittlesey S. N. , Research methods in biomechanics, Second edition. Champaign, Illinois: Human Kinetics, 2014. 51 Thelen D. G. and Anderson F. C. , “Using computed muscle control to generate forward dynamic simulations of human walking from experimental data,” J. Biomech., vol. 39 , no. 6 , pp. 1107–1115, Jan. 2006, doi: 10.1016/j.jbiomech.2005.02.010 16023125 52 Neptune R. R. , McGowan C. P. , and Kautz S. A. , “Forward Dynamics Simulations Provide Insight Into Muscle Mechanical Work During Human Locomotion,” Exerc. Sport Sci. Rev., vol. 37 , no. 4 , pp. 203–210, Oct. 2009, doi: 10.1097/JES.0b013e3181b7ea29 19955870 53 Faber H. , van Soest A. J. , and Kistemaker D. A. , “Inverse dynamics of mechanical multibody systems: An improved algorithm that ensures consistency between kinematics and external forces,” PLOS ONE, vol. 13 , no. 9 , p. e0204575, Sep. 2018, doi: 10.1371/journal.pone.0204575 30265727 54 Davis P. J. and Rabinowitz P. , Methods of Numerical Integration. Courier Corporation, 2007. 55 Hairer E. and Wanner G. , Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Berlin, Heidelberg: Springer Berlin Heidelberg, 1991. 56 Gavrea B. I. , Anitescu M. , and Potra F. A. , “Convergence of a Class of Semi-Implicit Time-Stepping Schemes for Nonsmooth Rigid Multibody Dynamics,” SIAM J. Optim., vol. 19 , no. 2 , pp. 969–1001, Jan. 2008, doi: 10.1137/060675745 57 Potra F. A. , “Runge-Kutta Integrators for Multibody Dynamics,” Mech. Struct. Mach., vol. 23 , no. 2 , pp. 181–197, Jan. 1995, doi: 10.1080/08905459508905234 58 (Edwin) Chen Y. J. , Sheen S. H. , Ascher U. M. , and Pai D. K. , “SIERE: A Hybrid Semi-Implicit Exponential Integrator for Efficiently Simulating Stiff Deformable Objects,” ACM Trans. Graph., vol. 40 , no. 1 , pp. 1–12, Feb. 2021, doi: 10.1145/3410527 59 Hairer E. , Roche M. , and Lubich C. , “Runge-Kutta methods for differential-algebraic equations,” in The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods, Hairer E. , Roche M. , and Lubich C. , Eds., in Lecture Notes in Mathematics. Berlin, Heidelberg: Springer, 1989, pp. 14–22. doi: 10.1007/BFb0093949 60 Pàmies-Vilà R. and Font-Llagunes J. , “Validation of the Inverse Dynamic Analysis of Human Gait Using a Forward Dynamics Approach,” presented at the Proceedings of the ASME Design Engineering Technical Conference, Aug. 2013. doi: 10.1115/DETC2013-13023 61 Jang S. J. and Choi Y. J. , “SYSTEM AND METHOD FOR EXECUTING A SIMULATION OF ACONSTRAINED MULTI-BODY SYSTEM,” p. 28, 2017. 62 Feldman A. , “Change in muscle length as a result of equilibrium changes in the muscle-load system,” Biofizika, vol. 19 , no. 3 , pp. 534–538, Jun. 1974.4422366 63 Bizzi E. , “Central and Peripheral Mechanisms in Motor Control,” in Advances in Psychology, Elsevier, 1980, pp. 131–143. doi: 10.1016/S0166-4115(08)61942-1 64 Hogan N. , “Impedance Control: An Approach to Manipulation,” in 1984 American Control Conference, San Diego, CA, USA: IEEE, Jul. 1984, pp. 304–313. doi: 10.23919/ACC.1984.4788393 65 Sobinov A. , “Description of motor control using inverse models,” PhD, West Virginia University Libraries, 2019. doi: 10.33915/etd.4073 66 Tucker M. R. et al. , “Control strategies for active lower extremity prosthetics and orthotics: a review,” J. NeuroEngineering Rehabil., vol. 12 , no. 1 , p. 1, Jan. 2015, doi: 10.1186/1743-0003-12-1 25557982 67 Buczek F. L. , Kepple T. M. , Siegel K. L. , and Stanhope S. J. , “Translational and rotational joint power terms in a six degree-of-freedom model of the normal ankle complex,” J. Biomech., vol. 27 , no. 12 , pp. 1447–1457, Dec. 1994, doi: 10.1016/0021-9290(94)90194-5 7806552