Efficient trajectory optimization for curved running using a 3D musculoskeletal model with implicit dynamics

Trajectory optimization with musculoskeletal models can be used to reconstruct measured movements and to predict changes in movements in response to environmental changes. It enables an exhaustive analysis of joint angles, joint moments, ground reaction forces, and muscle forces, among others. However, its application is still limited to simplified problems in two dimensional space or straight motions. The simulation of movements with directional changes, e.g. curved running, requires detailed three dimensional models which lead to a high-dimensional solution space. We extended a full-body three dimensional musculoskeletal model to be specialized for running with directional changes. Model dynamics were implemented implicitly and trajectory optimization problems were solved with direct collocation to enable efficient computation. Standing, straight running, and curved running were simulated starting from a random initial guess to confirm the capabilities of our model and approach: efficacy, tracking and predictive power. Altogether the simulations required 1 h 17 min and corresponded well to the reference data. The prediction of curved running using straight running as tracking data revealed the necessity of avoiding interpenetration of body segments. In summary, the proposed formulation is able to efficiently predict a new motion task while preserving dynamic consistency. Hence, labor-intensive and thus costly experimental studies could be replaced by simulations for movement analysis and virtual product design.


S2 System Dynamics
The system dynamics were described implicitly as a function f() of the states x, the state derivativesẋ, and the controls u (Eq. 1). The function f() contained • identitiesq − dq dt = 0 for each degree of freedom (DOF), • multibody dynamics for each DOF (Eq. S1), • activation dynamics for each muscle tendon unit (MTU) (Eq. S2), and • contraction dynamics for each MTU (Eq. S3).

S2.1 Multibody Dynamics
The multibody dynamics were defined as follows: where q contained the global position, the global orientation, and the joint angles,q contained the global velocities and joint angular velocities, M(q) was the mass matrix, C(q,q) contained the Coriolis forces, G(q) contained the gravity forces, and J c was the Jacobian of the contact forces F c . τ was the sum of active joint torques generated by the muscles τ mus (Eq. S12), passive joint torques τ pas (Eq. S14), and joint torques due to external actuation torques τ ext (Eq. S15). The talus was assumed to be weightless to save computation within the multibody dynamics. This assumption can be made since the talus is a small body which cannot move independently from the other segments since no MTU is connected to it.

S2.2 Muscle Dynamics
The model was operated using 92 MTUs (see Table S1). MTUs were modeled as three element Hill-type muscles with a contractile element (CE), a parallel elastic element (PEE), and a series elastic element (SEE). The muscle dynamics were described as a function of the activation state a and the CE length state s: where the state variable s = l CE cos(φ ) denoted the projection of CE length on the muscle line of action for a specific pennation angle φ 3 ,ȧ denoted the time derivative of the activation a, n e denoted the neural excitation, and r(n e ) denoted the activation dynamics, which were determined as follows: The activation time constant T act = 10 ms and the deactivation time constantand T deact = 40 ms were identically to Hamner's model 1 . In Eq. S3, F SEE denoted the force in the SEE, F CE denoted the force in the CE, and F PEE was the force in the PEE. The force in the CE was determined by where F ISO was the maximum isometric force, f FL (l CE ) denoted the force-length relationship, and f FV (v CE ) the force-velocity relationship. The force-length relationship was described as where l CE was the length of the CE normalized to the optimal fibre length l CE,opt , w a width parameter of the muscle, equal to the square root of the shape factor used by Thelen 4 . The force-velocity relationship was described as where v CE was the normalized CE velocity, λ = 0.5025 + 0.5341a 5 modeled the activation dependence of the normalized maximum shortening velocity v CE,max = 10 l CE,opt /s 4 . g max = 1.8 was the maximum force amplification during lengthening 4 , A = 0.25 was the Hill curve parameter 6 , and c FV was determined as follows 7 : The force in the SEE and PEE were modeled as non-linear springs 7 : where l denoted the length of the element normalized to optimal fiber length l CE,opt . It was equal to l MTU (q) − l CE for the SEE, and equal to s for the PEE. l slack was the slack length normalized to optimal fiber length l CE,opt , k 1 = 1 was a small linear stiffness, which was added to aid the optimization such that the derivative with respect to the model states was never zero, and k 2 was the non-linear stiffness, which was equal to the following for the PEE and SEE, respectively: where h = 0.04 was the strain of the muscle at isometric force 8, 9 .

S2.3 Muscle-Joint Coupling
The torque in DOF j generated by muscle i was determined using the following equation: where l MTU,i (q) denoted the normalized muscle-tendon length depending on the current pose defined by the generalized coordinates q, and ∂ l MTU,i (q) ∂ q j denoted the muscle moment arm. A constant muscle moment arm would not have been accurate enough in the 3D model 3 . Hence, a polynomial function was fitted to describe the muscle-tendon length depending on the joint angles to match the moment arms available in OpenSim, since polynomials have well defined derivatives. The order of the polynomial was chosen such that the root mean square error in the moment arms was less than 5 % of the maximum moment arm for each muscle. The maximum possible order was set to four. Additionally, the range of motion used to determine the polynomial was reduced if the muscles wrapped around the bones incorrectly for large joint angles. Linear interpolation was used outside of this range of motion.

S2.4 Passive Joint Torques
Passive torques were added to the joint torques when the joint angle was outside of the normal range of motion. These torques were determined as follows: where K 2 = 5000 N m rad −2 . The range of motions, defined by q j,min and q j,max for the trunk and legs are the same as those used for the muscle moment arms. For the arms, the following ranges were used: arm flexion τ pas, j = τ pas,out, j − K 1 (q j − q j,neutral ) − Bq j , (S14)

3/8
where the stiffness was K 1 = 1 N m rad −1 and the damping was equal to B = 1 N m s rad −1 . q j,neutral was the neutral position of DOF j defined as default values in the OpenSim model file runMaD.osim (see supplementary material at www.simtk.org).

S2.5 Arm Torques
The DOF j of a arm was directly actuated by the torque τ ext, j = m j 10 N m , with torque control m j . For numerical reasons, a scaling was performed to obtain states and controls of same magnitude.

S2.6 Penetration-Based Ground Contact Model
Eight contact points were used at each foot to describe the contact with the ground. Four contact points were located at the toe segment and four at the calcaneus segment. In the OpenSim model, their location relative to the respective segment origin was defined using marker objects (see runMaD.osim). The vertical ground reaction force (GRF) in each contact point c was determined based on the ground penetration d: The horizontal GRFs in x-and z-directions were determined using a continuous approximation of the Coulomb friction: where µ k = 1 was the kinetic friction coefficient,ṗ c,x andṗ c,z were the sliding velocities of the contact point, andṗ c,x,0 = p c,z,0 = 10 −2 m s −1 was a small velocity parameter that ensured that the force was differentiable around zero velocity.

S3 Simulations
The bounds of states x and controls u used in the three optimization examples are summarized in Table S2. For standing, smaller ranges were chosen compared to running to avoid the simulation terminating in local optima.

S4 Results
Pelvis translation and joint moments of straight and curved running are shown in Figs. S1 and S2. Table S2. Bounds used to simulate standing, straight running, and curved running. For straight running, the pelvis position at the first node was fixed to q pel_tx [0] = 0 and q pel_tz [0] = 0. For curved running, the pelvis position at the first node was fixed to q pel_tx [0] = −r and q pel_tz [0] = 0. ∆t denotes the duration between two collocation nodes.   Figure S2. Pelvis translation and joint moments of the curved running simulation. The degrees of freedom (DOFs) are named according to their definition in the model file runMaD.osim. Black, red, and blue solid lines indicate the simulated variables of the torso, the right side, and left side, respectively. Shaded areas show mean ± standard deviation (SD) of inverse dynamics (ID) of the measured gait cycles of curved running. The horizontal pelvis translation cannot be directly compared to the measured data since the global frames were not aligned but rotated around the vertical axis.