Does joint impedance improve dynamic leg simulations with explicit and implicit solvers?

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.

2 23 the impeded behavior has not only kinematic but also kinetic attributes that require interpretation. 24 Biomechanical models could allow medical specialists to observe these dynamic variables and 25 instantaneously diagnose mobility issues that may otherwise remain unnoticed. However, the real-26 time and subject-specific dynamic computations necessitate the optimization these simulations. In 27 this study, we explored the effects of intrinsic viscoelasticity, choice of numerical integration 28 method, and decrease in sampling frequency on the accuracy and stability of the simulation. The  The central nervous system (CNS) evolved to orchestrate muscular and skeletal actions to 43 produce complex motor behaviors (1). For instance, in locomotion, intrinsic and sensory 44 feedback signals are integrated with descending visuomotor commands to fine-tune limb 45 stepping and foot placement (2,3). To explain this intricate interplay between neural and . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint 3 46 mechanical systems, previous research has produced a variety of biomechanical models and 47 analyzed their interactions with the environment (4). For example, a biomechanical model of the 48 arm was shown to explain the directional tuning property of the primary motor cortex neurons, 49 describing why some cortical neurons are most active in particular movement directions (5). 50 Similar models of neuromechanical interactions are believed to exist within neural pathways in 51 order to overcome transmission delays and nonlinear limb dynamics (6). The use of embedded 52 dynamical computations for motor control is illustrated in Fig. 1. The desired motion 53 transformed through the inverse dynamics yields control signals to muscles (Fig. 1A). The 54 oversight over the execution of these motor commands can also be monitored with the use of  Gait abnormalities can be diagnosed and treated more effectively using forward dynamic 59 simulations of body motion (7-9). However, subject-specific pathomechanics may be required 60 for effective rehabilitation (10,11). Simple models may capture the dynamics of the limb in real- 61 time, but they may fail to capture the muscle-specific actions that require musculoskeletal to derive an objective approximation of muscle dynamics to accurately describe arm and hand 67 muscles in the context of real-time simulations with reduced computational load (21,22). This . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made In this study, we developed a biomechanical framework for testing viscoelastic joint 92 constraints and numerical solver types to reduce simulation errors at low sampling rates. We 93 posited that the transient kinematic perturbations due to the large integration time step could be 94 smoothed by adding viscoelastic joint impedance to increase simulation stability. Furthermore, 95 we hypothesized that the viscoelastic compensation depends on the type of a used numerical 96 solver. We tested these hypotheses in a human bipedal model simulating locomotion. We used 97 three numerical solvers with varied combinations of joint impedance to simulate the inverse and 98 forward dynamics of the swing locomotor phase. We determine optimal viscous and elastic pair 99 that stabilized the kinematic solutions and produced close-to-physiological kinetics. In addition, 100 we evaluated feasibility of our approach for real-time applications.

103
This study aimed to test the stabilization effects of viscoelastic compensation, evaluate the 104 performance of numerical solvers, and identify sampling frequencies that produce accurate and 105 stable simulations of the leg dynamics. This was accomplished with the following four steps to    The hip, knee, and ankle joint angles, as well as the standing foot rotation with respect to the ground, 133 were calculated from the Dataset 1 using the biomechanics analysis package Visual3D (C-Motion, 134 Inc.). The data were low-pass filtered using a 2 nd order Butterworth filter with the cutoff frequency    (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint 8 160 as a combination of a revolute and universal 2 DOF primitives to avoid the gimbal lock problem 161 (undesired alignment of two perpendicular axes). The foot was rigidly connected to the ground (no 162 slip condition for the leg in stance) close to metatarsophalangeal (MTP) joint (75% of foot length), 163 and MTP was allowed to rotate in 3 DOFs. This modeled an approximate location of the foot center 164 of pressure during the single stance phase (40). Since the analysis was focused on the leg in swing, 165 the motion of leg in stance was determined by its kinematics in both forward and inverse simulations.

166
This hybrid design prevented kinetic errors from propagating from the standing to the swinging leg.

167
While 17 DOFs were simulated, we focused our analysis on the seven DOFs of the leg in swing.

169
A built-in recursive algorithm (Simulink R2022b) was used to simulate the inverse dynamics of the 170 modeled body; the explicit 4th-order Runge-Kutta and explicit/implicit Euler methods were used to 171 simulate the forward dynamics. Euler explicit integration is the least accurate and yet the fastest.

172
Euler implicit formulation is better suited to keep the integration error bound to the tolerated error; (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  186 Passive joint impedance was described with the rotational stiffness and damping (k,b) of the spring-187 dampers used to stabilize multibody simulations (42). An increase in the impedance decreased 188 angular errors, but also increased torque errors. This, therefore, presented an optimization problem.

197
(2) Optimization: (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made initially set to 5% was adjusted with the Holm-Bonferroni correction, and the difference between 208 simulations with and without the viscoelastic contribution was tested with Wilcoxon rank sum test 209 (Fig. 3-4). The linear regression demonstrated the correlations between the execution times and 210 simulation rates (Fig. 5). All analyses were performed using the Statistics and Machine Learning

214
In this study, we used viscoelastic joint impedance to reduce numerical inaccuracies in simulations 215 of limb dynamics. The impedance increasing from proximal to distal joints significantly decreased 216 errors of explicit integration. The beneficial impedance stabilization was negligible when used with 217 implicit Euler method, which outperformed explicit Euler and Runge-Kutta methods. The added 218 impedance generated expected increase in torque errors, however, stabilized kinematics at low 219 sampling rates. Overall, the simulation errors were below 1% of peak-to-peak kinematics, and below 220 5% of peak-to-peak kinetics. All simulations were executed faster than in real-time. 222 We modified viscoelastic joint impedance to determine the optimal values that minimize both kinetic 223 and kinematic errors in simulations of leg movements. High impedance values increase torque errors 224 ( Fig. 2A) and decrease kinematic errors (Fig. 2B); low impedances have the opposite effect. We, 225 therefore, chose the optimal stiffness and damping sets as generating minimal kinetic and kinematic 226 errors at a given sampling rate. These sets were different for different DOFs and increased in value 227 from proximal to distal joints along the limb, i.e., the viscoelastic contribution was the lowest at the . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

A. OPTIMAL VISCOELASTIC IMPEDANCE
The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint 11 228 hip and the highest at the ankle. These values are shown with the blue outline in Fig. 2. 229 Supplementary tables (S1-S2) contain this information for all DOFs in the model.

230
The improvements in kinematic accuracy with additional impedance were not uniform across three  Table   249 S1-S2 for details). The viscoelastic contribution was therefore effective for forward simulations with 250 explicit methods and negligible for inverse simulations with fixed low sampling rates. In general, . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023.   simulations. However, small errors can be tolerated to increase integration accuracy with the 294 explicit methods especially when simulation frequency is low (Fig. 5). The implicit integration 295 method outperformed all explicit methods and required no additional compensation from the . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made impedance. Since explicit methods are standard in physics engines, the use of joint impedance 297 could improve their numerical stability; however, implicit methods deliver the best results. This 298 is supported by the comparative demonstration of handling object interactions with standard 299 physics engines and MuJoCo engine, which utilizes the implicit integration of velocity (26).

301
Accurate dynamic simulations of human limbs in real time require an efficient integration method.

302
Analytical solutions provide the best computational efficiency; however, they may not exist for 303 complex mechanical systems requiring fixed-or varied-step numerical methods to simulate body 304 movement (51). Fixed-step methods are typically used to avoid a sample mismatch between input 305 data acquisition and simulation rates; however, they suffer from the accumulating truncation error.

306
Variable-step methods are challenging to implement in real time because they require a predictive 307 extrapolation of real-time inputs and outputs in simulation. The integration methods of 308 approximating ODE can be explicit or implicit. Explicit schemes compute the future state of the 309 system (t+Δt) from its current state. However, they are often unstable when simulating 'stiff' 310 dynamical systems with slow and fast changing dynamical variables (52). For example, the knee 311 swiftly extends during the mid-swing and flexes slowly at the end of the swing of a human gait cycle; 312 the time of ground contact also introduces a sharp transition between swing and stance. Implicit 313 ("backward") schemes compute approximations using both the current and future states (41). This 314 requires solving for unknowns using a root-finding algorithm (for example, Newton's method). The 315 additional computational complexity within implicit methods is relatively less taxing than the 316 advantage of accuracy, as shown in Fig. 5. 317 We selected the simplest implicit integrator-the implicit Euler method-to solve the forward 318 dynamics of the human leg in swing. Other integrators often used for physical simulations include . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. implicit Euler over the explicit Euler and RK methods (Fig. 3-4). A modern state-of-the-art physical 332 engine MuJoCo uses RK, implicit Euler, and semi-implicit Euler modified to be fully implicit in the 333 velocity integration (26). While this method was not directly analyzed, our results provide insight 334 into its performance advantages.

335
As we demonstrate in Fig. 3-6, the solver selection is essential for optimizing a solution to a specific 336 forward dynamic problem. To model human movement, we recommend considering implicit 337 formulations or enforcing the explicit methods with viscoelastic impedance to balance simulation 338 speed, stability, and accuracy ( Fig.5-6, Fig. S2). 340 Choosing an optimal simulation frequency is a speed-accuracy tradeoff problem. We solved it by 341 computing inverse-forward dynamics for a range of simulation frequencies and assessing accuracy . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

B. SIMULATION FREQUENCY
The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint 16 342 and execution times. The inverse and forward simulation were performed sequentially, and the 343 numerical errors from both simulations accumulated, as in a similar study (57). The simulation 344 frequency needs to be minimized with acceptable level of errors for closed-loop control applications. 345 We chose the error thresholds (<1% kinematic and <5% kinetic) based on our previous evaluation 346 of errors in the approximation of musculoskeletal dynamics (21,22). We found forward and inverse 347 computations to be real-time accurate at frequencies as low as 200 Hz (Fig. 3 and 5, also see Fig.   348 S2). In simulation of object manipulation, the physics engine MuJoCo could maintain stability even 349 at lower simulation frequencies (62.5 Hz) due to many elegant methodological choices (29). This 350 suggests further benefits of using hybrid solvers; however, the confounding implementation 351 overhead maybe difficult to evaluate in software packages. 353 Passive mechanical impedance can stabilize physical simulation. Modeled using joint spring-354 dampers, they generate viscoelastic forces to oppose abrupt changes in position and velocity (58).

355
However, in inverse simulations, additional impedance may lead to an overestimation of joint forces. 356 We used this tradeoff to constraint our optimization and find the joint impedance that both reduced 357 kinematic errors and minimized force distortion. Our analyses showed that mechanical stabilization 358 did not improve the accuracy of implicitly computed kinematics, but, in general, improved the results 359 of explicit method. As a result, the hypothesized increase in simulation stability and the reduction in 360 transient kinematic perturbations were true for the explicit but not for the implicit integrators. It is 361 important to note, however, that the impedance-mediated stabilization is similar to kinematic 362 filtering with the correction to the middle of range of motion set by spring rest length. Here, we 363 applied the standard recommendation for filtering kinematics (37) and observed benefits of this . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint 17 364 approach. It is likely that in simulations that include abrupt mechanical interactions the 365 impedance-mediated stabilization may augment both implicit and explicit methods.

366
The concept of impedance-mediated stabilization in our study is applied to reducing numerical 367 errors, but it is also theoretically similar to the physiological mechanism. Muscle intrinsic properties 368 and co-contraction generate joint impedance to resist external perturbations. Moreover, the 369 impedance was hypothesized to be the controlled variable for movement and posture control (59-   Hz sampling frequency) and no viscoelastic contribution was optimal. The simulation was 401 kinematically (error <1%) and kinetically (error <5% peak-to-peak) accurate in real-time. 402 However, the viscoelastic impedance was essential to achieve accurate forward dynamic 403 solutions with explicit solvers.  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023.     . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. 626 Abbreviations: f-e-flexion-extension; ab-ad-abduction-adduction; int-ext-internal-external 627 rotation; ever-inv-eversion-inversion.

628
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2023. ; https://doi.org/10.1101/2023.02.09.527805 doi: bioRxiv preprint