Supplementary Information
Efﬁcient trajectory optimization for curved running
using a 3D musculoskeletal model with implicit
dynamics
Marlies Nitschke1,*, Eva Dorschky1, Dieter Heinrich2, Heiko Schlarb3, Bjoern M.
Eskoﬁer1, Anne D. Koelewijn1,4, and Antonie J. van den Bogert4
1Machine Learning and Data Analytics Lab, Department of Computer Science, Friedrich-Alexander-Universit¨at
Erlangen-N¨urnberg (FAU), Erlangen, Germany
2Department of Sport Science, University of Innsbruck, Innsbruck, Austria
3adidas AG, Herzogenaurach, Germany
4Department of Mechanical Engineering, Cleveland State University, Cleveland, USA
*marlies.nitschke@fau.de
S1 Model Adaptations
The proposed “running model for motions in all directions”, short “runMaD”, was adapted from a model proposed by Hamner
et al.1 in order to simulate running with directional changes. We changed the sequence of pelvis rotations from tilt (z), list (x),
rotation (y) to rotation (y), obliquity (x), tilt (z) (see Fig. 1 for the global coordinate system). The updated rotation sequence is
in agreement with clinical understanding2. Hence, results are interpretable in clinical analysis independently of the movement
direction. Furthermore, this sequence simpliﬁes the speciﬁcation of directional tasks in simulations, e.g. in the presented curved
running simulation.
Pronation and supination of the foot was enabled by unlocking the subtalar joint. The metatarsophalangeal (mtp) joint was
also unlocked for roll over of the foot. For both joints, a range of motion from −90◦ to 90◦ was allowed. In order to ﬁt the
recorded data of fast running, the upper limit of knee ﬂexion was increased from 120◦ to 160◦. Additionally, the range of motion
of the pronation/supination angle at the elbow was enlarged from [0◦,90◦] to [0◦,150◦] and the default pronation/supination
angle was set to 90◦ such that the palms were pointing towards the body for zero torque. The default elbow ﬂexion was set to
5◦ to be within the range of motion.
S2 System Dynamics
The system dynamics were described implicitly as a function f() of the states x, the state derivatives ˙x, and the controls u
(Eq. 1). The function f() contained
• identities ˙q− 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 deﬁned as follows:
M(q) ¨q + C(q, ˙q) ˙q + G(q) − JcT Fc − τ = 0,
(S1)
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 Jc
was the Jacobian of the contact forces Fc. τ 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:
˙a − r(ne) (ne − a) = 0,
(S2)
FSEE(lMTU(q),s) − (FCE(a,s, ˙s) + FPEE(s)) cos(φ(s)) = 0,
(S3)
where the state variable s = lCE cos(φ) denoted the projection of CE length on the muscle line of action for a speciﬁc pennation
angle φ 3, ˙a denoted the time derivative of the activation a, ne denoted the neural excitation, and r(ne) denoted the activation
dynamics, which were determined as follows:
r(ne) =
ne
Tact
+ 1 − ne
Tdeact
.
(S4)
The activation time constant Tact = 10ms and the deactivation time constantand Tdeact = 40ms were identically to Hamner’s
model1. In Eq. S3, FSEE denoted the force in the SEE, FCE denoted the force in the CE, and FPEE was the force in the PEE.
The force in the CE was determined by
FCE = a fFL(lCE) fFV(vCE) FISO ,
(S5)
where FISO was the maximum isometric force, fFL(lCE) denoted the force-length relationship, and fFV(vCE) the force-velocity
relationship. The force-length relationship was described as
fFL(lCE) = exp
 
−
lCE − 1
w
2!
,
(S6)
where lCE was the length of the CE normalized to the optimal ﬁbre length lCE,opt, w a width parameter of the muscle, equal to
the square root of the shape factor used by Thelen4. The force-velocity relationship was described as
fFV(vCE) =
(
λ vCE,max + vCE
λ vCE,max − vCE / A
if vCE < 0
gmax vCE + λ cFV
vCE + λ cFV
if vCE ≥ 0,
(S7)
where vCE was the normalized CE velocity, λ = 0.5025 + 0.5341a5 modeled the activation dependence of the normalized
maximum shortening velocity vCE,max = 10lCE,opt/s4. gmax = 1.8 was the maximum force ampliﬁcation during lengthening4,
A = 0.25 was the Hill curve parameter6, and cFV was determined as follows7:
cFV = vCE,max A (gmax − 1)
A + 1
.
(S8)
The force in the SEE and PEE were modeled as non-linear springs7:
F(l) =
(
FISO k1 (l − lslack)
if l ≤ lslack
FISO k1 (l − lslack) + FISO k2 (l − lslack)2
if l > lslack ,
(S9)
where l denoted the length of the element normalized to optimal ﬁber length lCE,opt. It was equal to lMTU(q) − lCE for the
SEE, and equal to s for the PEE. lslack was the slack length normalized to optimal ﬁber length lCE,opt, k1 = 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
k2 was the non-linear stiffness, which was equal to the following for the PEE and SEE, respectively:
k2,PEE =
1
w2 ,
(S10)
k2,SEE =
1
(h lslack,SEE)2 ,
(S11)
where h = 0.04 was the strain of the muscle at isometric force8,9.
2/8
Table S1. Abbreviations of muscles used in the model “runMaD”.
Abbreviation
Muscle Name
glut_min/med/max
gluteus minimus/medius/maximus
semimem
semimembranosus
semiten
semitendinosus
bifemsh/bifemlh
biceps femoris (short/long head)
sar
sartorius
add_long/brev/mag
adductor longus/brevis/magnus
tﬂ
tensor fascia latae
pect
pectineus
grac
gracilis
iliacus
iliacus
psoas
psoas
quad_fem
quadratus femoris
gem
gemellus
peri
piriformis
rect_fem
rectus femoris
vas_med/int/lat
vastus medialis/intermedius/lateralis
med/lat_gas
gastrocnemius (medial/lateral head)
soleus
soleus
tip_post/ant
tibialis posterior/anterior
ﬂex_dig/hal
ﬂexor digitorum/hallucis
per_long/brev/tert
peroneus longus/brevis/tertius
ext_dig/hal
extensor digitorum/hallucis
ercspn
erector spinae
intobl/extobl
internal/external abdominal oblique
S2.3 Muscle-Joint Coupling
The torque in DOF j generated by muscle i was determined using the following equation:
τmus,i,j = − ∂ lMTU,i(q)
∂ qj
FSEE,i (lMTU,i(q),lCE,i),
(S12)
where lMTU,i(q) denoted the normalized muscle-tendon length depending on the current pose deﬁned by the generalized
coordinates q, and ∂ lMTU,i(q)
∂ qj
denoted the muscle moment arm. A constant muscle moment arm would not have been accurate
enough in the 3D model3. Hence, a polynomial function was ﬁtted to describe the muscle-tendon length depending on the joint
angles to match the moment arms available in OpenSim, since polynomials have well deﬁned 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. The following ranges were used: hip ﬂexion [−28◦,78◦], hip adduction [−13◦,13◦], hip
rotation [−3◦,18◦], knee angle [−118◦,8◦], ankle angle [−38◦,38◦], subtalar angle [−13◦,13◦], mtp angle [−8◦,48◦], lumbar
extension [−38◦,3◦], lumbar bending [−8◦,8◦], and lumbar rotation [−18◦,18◦].
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:
τpas,out, j =
(
K2 (qj −qj,min)2
if q j < q j,min
−K2 (qj −qj,max)2
if q j > q j,max ,
(S13)
where K2 = 5000Nmrad−2. The range of motions, deﬁned by qj,min and qj,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 ﬂexion [−88◦,88◦], arm adduction
[−118◦,88◦], arm rotation [−88◦,88◦], elbow ﬂexion [2◦,148◦], and pronation/supination [2◦,148◦]. For numerical reasons, a
small stiffness was used for the full range of motion, such that the derivative of the joint moment with respect to the joint angle
was never zero:
τpas, j = τpas,out, j −K1 (qj − q j,neutral) − B ˙qj ,
(S14)
3/8
where the stiffness was K1 = 1Nmrad−1 and the damping was equal to B = 1Nmsrad−1. qj,neutral was the neutral position of
DOF j deﬁned as default values in the OpenSim model ﬁle 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 10Nm,
(S15)
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 deﬁned 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:
Fc,y(d) = k d (1 − b ˙pc,y),
(S16)
where k = 100BWm−1 was the stiffness of the ground normalized to body weight (BW) per meter. The damping constant
b = 0.75sm−110 was multiplied with the vertical velocity of the contact point ˙pc,y. The ground penetration d was determined
from the vertical position of the contact point pc,y and the size of the transition region pc,y,0 = 10−3 m between contact and no
contact:
d = 1
2
q
p2c,y + p2
c,y,0 − pc,y

.
(S17)
The horizontal GRFs in x- and z-directions were determined using a continuous approximation of the Coulomb friction:
Fc,x(Fc,y, ˙pc,x) = − µk Fc,y
˙pc,x
q
˙p2c,x + ˙p2
c,x,0
,
(S18)
Fc,z(Fc,y, ˙pc,z) = − µk Fc,y
˙pc,z
q
˙p2c,z + ˙p2
c,z,0
,
(S19)
where µk = 1 was the kinetic friction coefﬁcient, ˙pc,x and ˙pc,z were the sliding velocities of the contact point, and ˙pc,x,0 =
˙pc,z,0 = 10−2 ms−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.
4/8
Table S2. Bounds used to simulate standing, straight running, and curved running. For straight running, the pelvis position at
the ﬁrst node was ﬁxed to qpel_tx[0] = 0 and qpel_tz[0] = 0. For curved running, the pelvis position at the ﬁrst node was ﬁxed to
qpel_tx[0] = −r and qpel_tz[0] = 0. ∆t denotes the duration between two collocation nodes.
Parameter
Unit
Standing
Running
Lower
Upper
Lower
Upper
qpelvis_rotation
degree
0
0
-90
90
qpelvis_obliquity
degree
-5
5
-90
90
qpelvis_tilt
degree
-30
30
-90
90
qpel_tx
m
0
0
-5
7
qpel_ty
m
0.5
1.5
-1
2
qpel_tz
m
0
0
-3
3
qhip_flexion
degree
-30
30
-120
120
qhip_adduction
degree
-10
10
-120
120
qhip_rotation
degree
-30
30
-120
120
qknee_angle
degree
-30
10
-160
10
qankle_angle
degree
-30
30
-90
90
qsubtalar_angle
degree
-30
30
-90
90
qmtp_angle
degree
-30
30
-90
90
qlumbar_extension
degree
-10
10
-90
90
qlumbar_bending
degree
-10
10
-90
90
qlumbar_rotation
degree
-10
10
-90
90
qarm_flex
degree
-40
40
-40
40
qarm_add
degree
-40
40
-40
40
qarm_rot
degree
-40
40
-40
40
qelbow_flex
degree
0
150
0
150
qpro_sup
degree
0
150
0
150
˙q
rads−1 or ms−1
-30
30
-30
30
s
-
0
5
0
5
a
-
0
5
0
5
ne
-
0
5
0
5
m
-
-5
5
-5
5
5/8
0
100
0
3.5
Translation in m
pelvis_tx
0
100
0.95
1.15
pelvis_ty
0
100
−0.1
0.1
pelvis_tz
0
100
−250
150
Moment in Nm
hip_ﬂexion
0
100
−250
100
hip_adduction
0
100
−70
70
hip_rotation
0
100
−150
400
knee_angle
0
100
−300
100
Moment in Nm
ankle_angle
0
100
−40
120
subtalar_angle
0
100
−30
30
mtp_angle
0
100
−100
300
lumbar_extension
0
100
−250
150
Moment in Nm
lumbar_bending
0
100
−100
100
lumbar_rotation
0
100
−30
30
arm_ﬂex
0
100
−30
30
Gait Cycle in %
arm_add
0
100
−30
30
Gait Cycle in %
Moment in Nm
arm_rot
0
100
−30
30
Gait Cycle in %
elbow_ﬂex
0
100
−30
30
Gait Cycle in %
pro_sup
Simulated
Measured: Mean ± SD
Torso
Right side
Left side
Figure S1. Pelvis translation and joint moments of the straight running simulation. The degrees of freedom (DOFs) are
named according to their deﬁnition in the model ﬁle 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 straight running.
6/8
0
100
−4
3
Translation in m
pelvis_tx
0
100
0.9
1.1
pelvis_ty
0
100
0
2
pelvis_tz
0
100
−250
150
Moment in Nm
hip_ﬂexion
0
100
−250
100
hip_adduction
0
100
−70
70
hip_rotation
0
100
−150
400
knee_angle
0
100
−300
100
Moment in Nm
ankle_angle
0
100
−40
120
subtalar_angle
0
100
−30
30
mtp_angle
0
100
−100
300
lumbar_extension
0
100
−250
150
Moment in Nm
lumbar_bending
0
100
−100
100
lumbar_rotation
0
100
−30
30
arm_ﬂex
0
100
−30
30
Gait Cycle in %
arm_add
0
100
−30
30
Gait Cycle in %
Moment in Nm
arm_rot
0
100
−30
30
Gait Cycle in %
elbow_ﬂex
0
100
−30
30
Gait Cycle in %
pro_sup
Simulated
Measured: Mean ± SD
Torso
Right side
Left side
Figure S2. Pelvis translation and joint moments of the curved running simulation. The degrees of freedom (DOFs) are named
according to their deﬁnition in the model ﬁle 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.
7/8
References
1. Hamner, S., Seth, A. & Delp, S. L. Muscle contributions to propulsion and support during running. J. Biomech. 43,
2709–2716 (2010).
2. Baker, R. Pelvic angles: a mathematically rigorous deﬁnition which is consistent with a conventional clinical understanding
of the terms. Gait & Posture 13, 1–6 (2001).
3. Van den Bogert, A. J., Blana, D. & Heinrich, D. Implicit methods for efﬁcient musculoskeletal simulation and optimal
control. Procedia IUTAM 2, 297–316 (2011).
4. Thelen, D. G. Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. J.
Biomech. Eng. 125, 70 (2003).
5. Chow, J. W. & Darling, W. G. The maximum shortening velocity of muscle should be scaled with activation. J. Appl.
Physiol. 86, 1025–1031 (1999).
6. Winters, J. M. An improved muscle-reﬂex actuator for use in large-scale neuromusculoskeletal models. Annals Biomed.
Eng. 23, 359–374 (1995).
7. McLean, S. G., Su, A. & Van den Bogert, A. J. Development and validation of a 3-D model to predict knee joint loading
during dynamic movement. J. biomechanical engineering 125, 864–874 (2003).
8. Van Soest, A. J. & Bobbert, M. F. The contribution of muscle properties in the control of explosive movements. Biol.
cybernetics 69, 195–204 (1993).
9. Miller, R. H. Hill-based muscle modeling. In Handbook of Human Motion, 373–394 (Springer, 2018).
10. Gerritsen, K. G., van den Bogert, A. J. & Nigg, B. M. Direct dynamics simulation of the impact phase in heel-toe running.
J. Biomech. 28, 661–668 (1995).
8/8

Efficient trajectory optimization for curved running using a 3D musculoskeletal model with implicit dynamics.
Nitschke, Marlies,Dorschky, Eva,Heinrich, Dieter,Schlarb, Heiko,Eskofier, Bjoern M,Koelewijn, Anne D,van den Bogert, Antonie J