1
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  
| https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports
Efficient trajectory optimization 
for curved running using 
a 3D musculoskeletal model 
with implicit dynamics
Marlies Nitschke1*, Eva Dorschky1, Dieter Heinrich2, Heiko Schlarb3, Bjoern M. Eskofier1, 
Anne D. Koelewijn1,4 & Antonie J. van den Bogert4
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.
In recent years, interest in musculoskeletal simulation to reconstruct and predict human movements has been 
 growing1. Motion reconstruction based on captured data yields insight into further variables of interest, e.g. joint 
moments or muscle  forces2–7. Furthermore, simulations can be applied to predict changes of kinematics as well 
as joint and muscle function in response to interventions or environmental changes. They can support decisions 
in orthopaedic  surgeries8,9 and the design of  prostheses10,11,  exoskeletons12, or  shoes13. Therefore, predictive 
simulations can replace time-consuming and expensive prototyping and experimental studies.
Commonly, biomechanical parameters are computed in a consecutive approach using inverse kinematics (IK), 
inverse dynamics (ID), and static optimization. This results in inconsistencies between kinematics and  kinetics14. 
Additionally, each time step is analyzed separately in a discrete set of optimization problems rather than solving 
one optimization problem over time. Since either the motion or the forces are not simulated but prescribed, 
these methods cannot be applied to predict novel movements. These limitations can be overcome by solving 
an open-loop optimal control problem, also known as trajectory optimization. Movements are reconstructed 
or predicted by obtaining state and control trajectories from one single constrained non-linear optimization 
 problem7. Alternatively, human motion can be simulated using a neuromuscular model with a controller based 
on  reflexes15,16 or central pattern  generators17.
Because of their relevance to daily activities and sports, it is important to also simulate movements with 
changes in direction additional to straight walking and running. It is also important that such movements can 
be optimized with respect to performance and/or injury risk, without having to collect human motion data from 
those movements. This will expand the use of existing data and avoid potentially risky human experiments. 
OPEN
1Machine Learning and Data Analytics Lab, Department of Computer Science, Friedrich-Alexander-Universität 
Erlangen-Nürnberg (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. *email: marlies.nitschke@fau.de
2
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Potential applications are in knee injury  prevention18, shoe performance in curved  running19 or 3D controllers 
for exoskeletons and active  prostheses20. Simulations of gait including turning were simulated with a reflex-
based  controller16, but reflex loops limit the space of possible control inputs. Trajectory optimizations find 
open-loop control trajectories without limiting the inputs, but have only been performed in two dimensional 
(2D)  space2,7,11,13,21–23 or were restricted to forward  motions24–27. A three dimensional (3D) model is needed 
to simulate movements with directional changes leading to a high-dimensional solution space of the optimal 
control problem. Therefore, effective numerical methods are necessary to put 3D optimal control simulations 
into application. Specifically, we need the capability to find solutions from an initial guess that is far from the 
solution, the solution should be found in a reasonable amount of time, and the solutions should exactly satisfy 
task requirements, such as a specified running speed and change of direction. The latter is of critical importance 
for sports applications.
Often, trajectory optimization problems were solved using an explicit formulation of the multibody dynamics 
leading to large computational  costs25–30. The reconstruction of one cycle of walking or running took for example 
around 2 h for a 3D model with 21 degrees of freedom (DOFs) and 66 muscle tendon units (MTUs)25 using 
direct collocation with OpenSim’s explicit implementation of the  dynamics5,31. This seems to be inefficient when 
taking into account that the initial guess was close to the simulation results since it was generated from human 
motion data of the same movement  task25.
An implicit formulation of the model dynamics was proposed by Van den Bogert et al.7 for a 2D model to 
improve the numerical conditioning of the optimal control problem and thus reduced computational cost. They 
used direct collocation to solve tracking as well as predictive optimal control simulations. This approach was 
successfully applied to analyze loading asymmetry in transtibial amputee  gait11 and the effect of midsole materi-
als of  shoes13. Recently, Falisse et al.24 used implicit dynamics to develop a framework for rapid simulations of 
a 3D model. They generated predictive simulations of walking and running with a full-body 3D model with 29 
DOFs and 92  MTUs32 on average in 36 min. However, until now no movements with directional changes in 3D 
space were simulated.
The purpose of our work was to further extend the current state-of-the-art in trajectory optimization for 
musculoskeletal models by computational efficient simulations of movements with directional changes. To this 
end, we created a complex full-body 3D musculoskeletal model called “running model for motions in all direc-
tions”, short “runMaD”, adapted from Hamner et al.32. To reduce computational cost, dynamics were formulated 
implicitly, and derivatives were formulated analytically. We demonstrated the efficacy, the tracking capabilities, 
as well as the predictive power of the proposed trajectory optimization with implicitly formulated dynamics 
and direct collocation using three simulations: prediction of static standing, tracking of straight running, and 
prediction of curved running.
Methods
In the following, we describe the developed musculoskeletal model, the general trajectory optimization approach 
and how we generated simulations of standing, straight running, and curved running.
3D musculoskeletal model. 
The proposed musculoskeletal model “runMaD” is a full-body 3D model 
with 33 DOFs, operated using 92 MTUs in the trunk and legs and 10 torque actuators in the arms (see Fig. 1). 
This model was adapted from the OpenSim model created by Hamner et al.32. The order of rotations in the pelvis 
was  changed33 and the subtalar and metatarsophalangeal (mtp) joints were unlocked to simulate movements 
with directional changes. Ranges of motion were enlarged for knee flexion and pronation/supination angle at the 
elbow to fit the recorded motion. Muscular and segmental properties were taken from Hamner’s model. We refer 
to section S1 in the Supplementary Information for a detailed description of all model adaptations.
All muscles were modeled as three element Hill-type muscles with a contractile element (CE) with contrac-
tion and activation dynamics, a parallel elastic element (PEE), and a series elastic element (SEE) (Fig. 2)34. The 
dynamics were described implicitly for each muscle with respect to the activation a and the state variable s , which 
was the projection of the CE length lCE on the muscle line of  action7. The variable s was used instead of lCE to 
avoid singularities with respect to the pennation angle φ7. Muscle-tendon lengths were described as polynomial 
functions of joint angles, which were fitted using the muscle moment arm data of Hamner et al.32. In accordance 
with Falisse et al.24, polynomial functions were chosen since they have well-defined derivatives.
The model’s state vector was defined by x = ( q
˙q
s
a )T , where q contained the DOFs, ˙q the derivatives 
of the DOFs, s the CE length state of all muscles, and a the activation state of all muscles. The control vector 
was defined as u = ( ne
m )T , where ne denoted the neural excitation of all muscles and m the arm actuation 
torque divided by 10 Nm for each of the DOFs in the arm (Eq. S15 in Supplementary Information). Dynamics 
of the musculoskeletal system were combined with a penetration-based ground contact model to describe the 
full dynamics of the model implicitly as function of the states x , the state derivatives ˙x , and the controls u:
Details on the system dynamics are given in section S2 in the Supplementary Information.
Trajectory optimization. 
Optimal control problems were formulated to generate movement simulations. 
The goal was to find a state trajectory x(t) , a control trajectory u(t) , and the duration of the simulated movement 
Tsim such that the objective function J(x(t), u(t), Tsim) was minimized with respect to the following constraints:
(1)
f(x(t), ˙x(t), u(t)) = 0.
(2)
f(x(t), ˙x(t), u(t)) = 0
(dynamic equilibrium)
3
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
The periodicity constraint ensured that the state of the model at the end of the gait cycle, Tsim , was equal to the 
state at the beginning rotated in the horizontal plane with Rper and shifted by horizontal translation tper.
The objective was defined as weighted sum of a tracking term, a muscular effort term, a torque term, and a 
regularization term:
(3)
xL(t) ≤ x(t) ≤ xU(t)
(bounds on states)
(4)
uL(t) ≤ u(t) ≤ uU(t)
(bounds on controls)
(5)
x(Tsim) = Rper x(0) + tper
(periodicity of states)
(6)
u(Tsim) = u(0)
(periodicity of controls)
(7)
J(x(t), u(t), Tsim) = WTrack Jtrack + Wmus Jmus + Wtor Jtor + Wreg Jreg.
acromial
elbow
ulna
radius
hand
torso
pelvis
wrist
radioulnar
humerus
femur
tibia
calcaneus
talus
toes
hip
knee
ankle
subtalar
metatarso-
phalangeal
y
G
x
G
Figure 1.  Musculoskeletal model “runMaD” with segments in black, joints in blue, ground contact points in 
pink, and the global coordinate system in green. The musculoskeletal model was visualized using OpenSim 4.0 
(https ://opens im.stanf ord.edu).
PEE
SEE
s
CE
ϕ
lCE
lMTU
Figure 2.  Hill-type muscle  model34 with contractile element (CE), parallel elastic element (PEE), series elastic 
element (SEE), muscle-tendon length lMTU , length of the CE lCE , and pennation angle φ . The state variable s 
represents the projection of lCE on the muscle line of action.
4
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
For tracking, the squared difference between simulated data ysim and the corresponding mean measured data 
µymeas of multiple gait cycles was minimized for all time points t:
The squared difference was normalized by the variance of measured data σ 2
ymeas to make it dimensionless. The 
variance σ 2
ymeas was adapted to be at least 10 % of the mean of the variance to avoid division by small numbers. This 
is for example necessary for the variance of the ground reaction force (GRF) during swing phase. Furthermore, 
the terms were normalized by the duration of the simulation Tsim and weighted with WVar,i for each signal vari-
able i. Volume-weighted and cubed neural excitation ne was minimized for each of the Nmus muscles to reduce 
muscular effort and to solve the muscle redundancy problem:
where the ratio of muscle volume was computed with the maximum isometric force FISO and the optimal length 
of the CE lCE,opt . The muscle effort was divided by the cubic norm of horizontal translation speeds 
 vx
vz
3 to 
compensate for different running speeds. It was shown previously that muscle activation is linear to movement 
 speed35,36. Besides the control of the muscles, the torque controls mi actuating the arms were minimized:
Finally, a small regularization term was added to enhance convergence by minimizing the derivatives of the 
states x and controls u:
We used regularization for both states and controls since we found that this yields the lowest number of iterations 
without losing simulation accuracy.
Simulations. 
We performed three simulations: prediction of static standing, tracking of straight running, 
and prediction of curved running. In the following, we describe the data acquisition, the optimal control prob-
lems of the three simulations, and the details of implementation and solution process. An overview of the pipe-
line is given in Fig. 3.
Experimental data. We recorded straight and curved running of a male subject (92 kg, 1.95 m) using 42 reflec-
tive markers, 16 infrared cameras (Vicon MX, Oxford, UK), and two force plates (Kistler Instruments Corp, 
Winterthur, CH) for tracking and as reference. The sampling frequency was set to 200 Hz and 1000 Hz, respec-
tively. Straight running was performed at a speed of vx = 4.0 ms−1 and vz = 0 . Curved running was performed 
in a circle with radius r = 3.7 m at a norm horizontal speed of 
 vx
vz
 = 2.7 ms−1 . The API of OpenSim 4.05 
was called within MATLAB (Mathwork, Natick, MA, USA) to scale the generic model using marker trajecto-
ries in neutral pose with arms besides the body (N-pose) and to compute IK and ID. For ID, joint angles were 
filtered within OpenSim with a 3rd order dual-pass low-pass Butterworth filter with a cut-off frequency set to 
15  Hz37. The GRFs were filtered with the same filter to avoid artifacts in the computed joint  moments38. After 
processing, single gait cycles were extracted from right to right heel strike using the minimum of the right heel 
marker. The mean and standard deviation (SD) of 12 gait cycles were computed for straight and curved running 
after linearly interpolating to the number of samples of the shortest cycles. The subject gave informed consent 
prior to participation. The study was approved by the ethical committee of the Friedrich-Alexander-Universität 
Erlangen-Nürnberg (Re.-No. 106_13 B). All methods were carried out in accordance with relevant guidelines 
and regulations.
Standing. The goal of the prediction of static standing was to find a neutral pose of the model in equilibrium 
without data tracking. As this simulation was independent of time, ˙x was set to 0 and no periodicity constraints 
(Eqs. 5 and 6) were used. The weights of the objective terms were chosen empirically for all simulations. For 
standing, Wmus = Wtor = 1 was set. All bounds of x and u are provided in Table S2 in the Supplementary Infor-
mation. The simulation was solved 50 times with different random initial guesses for states x and controls u 
to reduce the likelihood of ending up in a local minimum. The result with the lowest objective was chosen as 
solution.
(8)
Jtrack =
1
Tsim
Tsim

0
Ntrack

i=1
WVar,i
ysim,i(t) − µymeas,i(t)
σymeas,i(t)
2
dt.
(9)
Jmus =
1
Tsim Nmus
 vx
vz
3
Tsim

0
Nmus

i=1





FISO,i lCE,opt,i
Nmus

i=1
FISO,i lCE,opt,i
ne,i(t)3




 dt,
(10)
Jtor =
1
Tsim Ntor
Tsim

0
Ntor

i=1
mi(t)2 dt.
(11)
Jreg =
1
Tsim (Nstates + Ncontrols)
Tsim

0
Nstates

i=1
˙xi(t)2 +
Ncontrols

i=1
˙ui(t)2

dt.
5
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Straight running. Straight running was reconstructed by tracking straight running data. All joint angles, the 
global orientation of the pelvis and the GRFs of both feet in all directions were tracked (Eq. 8) similar to our 
approach in  2D7. A weighted arithmetic mean was used to balance the influence between joint angle tracking 
and GRF tracking independently of the number of signals: WVar,i =
Wang/GRF
NAng Wang + NGRF WGRF with Wang = 1 and 
WGRF = 5 and the numbers of tracked signals NAng and NGRF . The running speed used for weighting in the 
muscular effort term (Eq. 9) was computed from the tracking data. The weights of the objective terms were set 
to WTrack = 1 , Wmus = 103 , and Wtor = 1 such that after optimization the weighted objectives of tracking and 
effort were of same scale. The weight of the regularization term was small.
In the straight running simulation, the periodicity of the gait cycle was ensured by allowing only translation 
in the horizontal plane with prescribed running speeds (Eq. 5). Hence, Rper in Eq. (5) was the identity matrix 
and only a translation for the global pelvis position was applied:
Additionally, the states x and controls u were limited by lower and upper bounds (Eqs. 3 and 4). To define a 
global start position of the motion, the pelvis position at the first node was fixed to qpel_tx[0] = qpel_tz[0] = 0 . 
The standing solution was used as initial guess.
Curved running. Curved running was predicted by tracking straight running data and constraining the model 
to run in a circle. Only joint angles and vertical GRFs were tracked to allow a circular motion (Eq. 8). The norm 
horizontal speed was obtained from the reference data of curved running to weight the muscular effort term 
(Eq. 9). In contrast to straight running, Wmus and Wtor were increased by factor 10 to allow more deviation from 
the tracking data.
With help of the periodicity constraint (Eq. 5), we ensured that the model ran counterclockwise in a circle 
around the y-axis so that the left leg was on the inside. The circle was centered at 
 qpel_tx
qpel_tz

=
 0
0

 with central 
angle θ . All entries in the state vector x were constrained to be equal for t = 0 and t = Tsim except for the global 
pelvis position and rotation:
The norm horizontal speed and the radius of the measured curved running were used to obtain the central angle:
To define a global start position of the motion, the pelvis position at the first node was fixed to qpel_tx[0] = −r 
and qpel_tz[0] = 0 . The straight running solution was used as initial guess.
(12)


qpel_tx(Tsim)
qpel_ty(Tsim)
qpel_tz(Tsim)

 =


qpel_tx(0)
qpel_ty(0)
qpel_tz(0)

 +
 vx Tsim
0
vz Tsim

.
(13)



qpel_tx(Tsim)
qpel_ty(Tsim)
qpel_tz(Tsim)
qpel_rot(Tsim)


 =



cos(θ) 0 sin(θ) 0
0
1
0
0
−sin(θ) 0 cos(θ) 0
0
0
0
1






qpel_tx(0)
qpel_ty(0)
qpel_tz(0)
qpel_rot(0)


 +



0
0
0
θ


 .
(14)
θ = 2 arcsin
 vx
vz
 Tsim
2 r

.
Figure 3.  Processing pipeline. Work presented in this paper is highlighted in blue and OpenSim applications 
in red. The unscaled and scaled model, the experimental data, as well as the simulation results are provided in 
OpenSim file formats in the electronic supplementary material and at https ://simtk .org/proje cts/runma d.
6
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Implementation and solution process. The implicit multibody dynamics of the skeletal model were derived 
using Autolev (Symbolic Dynamics Inc., Sunnyvale, CA, USA) including Jacobian matrices generated by sym-
bolic differentiation. Muscle dynamics and the ground contact model were implemented in C. All dynamic 
equations were compiled as MEX-functions in MATLAB. The objectives and the task constraints as well as their 
analytic derivatives were coded in MATLAB.
The optimal control problems were solved for the scaled musculoskeletal model using direct collocation and 
backward Euler discretization. One collocation node was used for the static standing simulation. 50 collocation 
nodes were chosen for the running simulations since we found in preliminary simulations that 50, 100 and 200 
collocation nodes yielded similar results. Running data was linearly interpolated to 50 samples. The non-linear 
optimization problems were solved with IPOPT 3.12.339. Settings were adapted to terminate the optimization at 
a tolerance of  10−5, a constraint violation tolerance of  10−3, and a complementary tolerance of  10−3. All optimiza-
tions were run on one core of a workstation with a 3.2 GHz Xeon E5-1660v4 processor.
Results
In total, 50 predictive simulations of static standing using different random initial guesses, one tracking simula-
tion of straight running, and one predictive simulation of curved running were solved. CPU times and iterations 
required for optimization are summarized in Table 1. The entire process from generating standing from random 
initial guesses to the prediction of curved running took 1 h 17 min 28 s. The tracking of straight running using 
standing as initial guess required more iterations than the prediction of curved running using straight running 
as initial guess since standing is quite a different task than running.
The results corresponded to upright standing and natural running motion (Fig. 4 and videos in the electronic 
supplementary material). Joint angles and GRFs of straight running were close to the reference data of straight 
running, since the difference was generally less than one SD (Fig. 5). However, simulated knee flexion was smaller 
during the swing phase compared to the reference data. The difference between the right and left subtalar angle 
was well represented in the simulation of straight running even though subtalar and mtp angles deviated from 
the reference data more than one SD. The differences in GRFs were larger than one SD between 20% and 40% of 
the gait cycle. The estimated muscle activation patterns of the 18 largest muscles were similar to the electromyo-
graphy (EMG) measurements reported by Cappellini et al.40 for straight running (Fig. 5).
Predicted joint angles and GRFs of curved running matched the reference data of curved running but were not 
as similar as for the tracking simulation (Fig. 6). In particular, the ranges of motion were underestimated in the 
hip and knee flexion, knee flexion was smaller during stance, and maximum vertical GRFs were overestimated. 
The pelvis rotation and the horizontal GRFs cannot be compared directly for curved running since the global 
frames of simulation and reference were not aligned but rotated around the vertical axis. Muscle activations for 
curved running were similar to straight running but of less amplitude and less symmetric (Fig. 6). The model 
Table 1.  Solver performance for the simulations of standing, straight running, and curved running. For 
standing, the CPU times and number of iterations are listed for the chosen result which had the lowest 
objective of the 50 simulations. Additionally, the mean and sum for all 50 standing simulations is given.
CPU time in mm:ss
Iterations
Standing
Result
00:10
355
Mean
00:14
513
Sum
11:26
25,669
Straight running
46:22
1,065
Curved running
19:40
715
Figure 4.  Stick figures showing the simulation process. From left to right: Initial guess used for standing, result 
of standing, result of straight running, and result of curved running. For visualization, every fifth and every 25th 
node was plotted for straight and curved running, respectively, and the curved running result was extended to 
fill a whole circle.
7
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Figure 5.  Joint angles, GRFs, and muscle activations of the 18 largest muscles of the straight running simulation 
at 4.0  ms−1. GRFs are scaled to body weight (BW) and muscle activations are normalized to the peak activation 
of straight running. The degrees of freedom (DOFs) and muscles are named according to their definition in the 
model file runMaD.osim. An overview of all muscles is provided in Table S1 in the Supplementary Information. 
Black, red, and blue solid lines indicate the simulated variables of the torso, the right side, and left side, 
respectively. For GRFs and joint angles, shaded areas show mean ± standard deviation (SD) of the measured 
gait cycles of straight running. For muscles, shaded red areas show mean electromyography (EMG) data of the 
right side for running at 3.3  ms−1 reported by Cappellini et al.40 normalized to the maximum of the simulated 
activations of each muscle.
8
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
and data of the simulated and measured movements is provided in OpenSim file formats in the electronic sup-
plementary material and at https ://simtk .org/proje cts/runma d.
Discussion
In this work, we presented a 3D full-body musculoskeletal model adapted to running with directional changes 
and an implicit formulation of its dynamics for efficient trajectory optimization using direct collocation. We gen-
erated a predictive simulation of standing, a tracking simulation of straight running, and a predictive simulation 
with directional change, more precisely curved running, to demonstrate the efficacy, the tracking capabilities, 
and the predictive power of the approach.
Computational efficiencies must be compared with caution due to influences of computational power, imple-
mentation, the initial guess, and the choice of the numerical solution method, e.g. direct  collocation7,25,26, multiple 
 shooting41,42, and simulated  annealing27. Nevertheless, it can be concluded that the CPU times of approximately 
46 min and 20 min for straight and curved running, respectively, were small compared to simulations with 
explicit  dynamics25,26. We were able to solve tracking simulations of straight running faster than Lin and  Pandy25 
and Lin et al.26, while starting from an initial guess that was constructed without IK, static optimization, or 
computed muscle control (CMC). Furthermore, our model “runMaD” had a higher complexity resulting in 
approximately 20 % more unknowns, i.e. states and controls, per collocation node. The CPU times were similar 
to Falisse et al.24, who simulated straight running, confirming the advantages of implicit formulation of model 
dynamics. Falisse et al.24 used algorithmic derivatives whereas we implemented them analytically. Future studies 
are needed to compare the efficiency of algorithmic and analytical derivatives in optimal control simulations. 
Previously, it was shown that using algorithmic derivatives is faster than finite differences which are used by 
 OpenSim43. Machine learning approaches were recently investigated in the field of computer graphics to speed up 
musculoskeletal  simulation44,45. Jiang et al.44 learned a mapping from muscle-actuation space to joint-actuation 
space which would have to be retrained if model parameters are changing. Lee et al.45 used a two-stage deep 
reinforcement learning approach to simulate a full-body 3D musculoskeletal model. However, muscle activation 
was not included in the reward of the trajectory mimicking which might lead to non-optimal use of the muscles.
In all three simulations, kinematics were natural and matched the reference data. However, smaller knee 
flexion during the swing phase resulted in a lower foot clearance in both running simulations (Figs. 5 and 6). This 
movement pattern might be more energy efficient since increased foot clearance requires more effort. Subtalar 
and mtp angles deviated from the reference data and seemed to compensate each other’s error (Figs. 5 and 6). The 
reference data might be erroneous since the subject wore shoes and the foot segments are small which increases 
sensitivity to errors in marker positions and soft tissue artifacts. Additionally, simplifications in the foot model 
likely caused inaccuracies in the simulation. Although our proposed model is allowing subtalar and mtp motion 
in contrast to the model of Hamner et al.32, our foot model does also not reflect the fine foot structures. For a 
detailed analysis of foot motion, a finite element model of the foot, like developed by Akrami et al.46, could be 
incorporated into the musculoskeletal model. Instead of data tracking, movements could be reconstructed by 
constraining the movement  path25,26,42. Nevertheless, the bounds cannot be exceeded and are difficult to choose 
especially if a change in motion should be predicted.
In the predictive simulation of curved running, interpenetration of the legs occurred to reduce muscular effort 
(see video in electronic supplementary material). As a result of this, kinematics deviated slightly from reference 
data, i.e. ranges of motion were underestimated in the hip and knee flexion (Fig. 6). Interpenetration could be 
avoided by prescribing a minimal distance between joint  origins24 which would however require prior knowledge 
about the motion path. Alternatively, a constraint or an error term could be added to prevent intersection of 
segments similar to what was done in computer  graphics47. In both cases, the use of a stochastic  environment48 
would be beneficial to avoid segments moving close to each other. Similarly, the predictive simulation was lacking 
knee flexion during stance since it does not account for uncertainties while minimizing  effort24,49.
Even though GRFs of straight running were tracked well, impact of the initial contact was not distinct due 
to a fast progression towards the forefoot (Fig. 5). This was caused by the relatively coarse sampling with 50 
nodes. Additionally, the simple foot and contact model might have influenced especially the impact. For curved 
running, maximum vertical GRFs were overestimated compared to the reference data (Fig. 6) probably due to 
the interpenetration of the legs. Nevertheless, the outer leg, i.e. the right leg, showed a higher maximum vertical 
GRF compared to the inner leg in agreement with the reference data.
Joint moments were smoother in the simulation compared to the ID since the effort term has a smoothing 
effect (Figs. S1 and S2 in the Supplementary Information). ID depends on filtering of joint kinematics and  GRFs38 
whereas trajectory optimization benefits from a physics-based filtering. The mtp moments cannot be compared 
to the reference until late stance phase since GRFs were applied to the calcaneus in ID whereas we simulated 
ground contact at the calcaneus and toe segment. However, the mismatch of simulated and measured subtalar 
and mtp angles, contributed to the deviation in joint moments.
In contrast to static optimization, muscle activations and controls were computed while accounting for 
muscle and tendon dynamics. This is especially important for the analysis of fast movements, like  running50. 
The simulated muscle activations for straight running followed mainly the patterns of EMG  measurements40 
(Fig. 5) despite a higher velocity (4.0  ms−1 vs. approximately 3.3  ms−1). In comparison to straight running, peak 
muscle activations were smaller for curved running since running velocity was smaller (4.0  ms−1 vs. 2.7  ms−1) 
(Fig. 6). Furthermore, muscle activations were less symmetric due to the asymmetric movement task. We were 
able to predict curved running at 2.7  ms−1 while using tracking data of straight running at a different velocity 
of 4.0  ms−1. This confirms the predictive power of the trajectory optimization, because the target and tracking 
velocity did not have to match.
9
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Figure 6.  Joint angles, GRFs, and muscle activations of the 18 largest muscles of the curved running simulation 
at 2.7  ms−1. GRFs are scaled to body weight (BW) and muscle activations are normalized to the peak activation 
of straight running. The degrees of freedom (DOFs) and muscles are named according to their definition in the 
model file runMaD.osim. An overview of all muscles is provided in Table S1 in the Supplementary Information. 
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 the measured gait cycles of curved running.
10
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
Three general limitations of trajectory optimization have to be mentioned. First, it cannot be ensured that the 
global optimum of the optimal control problem was found. In the standing simulation, we used multiple initial 
guesses to minimize the risk of ending in a local optimum. Second, reported kinematics and kinetics were affected 
by the objective function and thus by the choice of weights of the objective terms. GRFs were weighted more in 
comparison to joint angles since the GRFs were the only tracked signals containing information about the forces 
within the body. For higher tracking weight WTrack in the straight running simulation, GRFs and joint angles were 
tracked better. However, this came at the cost of non-smooth activation signals. The signals contained alternating 
phases of activation and deactivation to allow the joint angles to closely match the data. When the same weights 
of the straight running simulation were used for the curved running simulation, i.e. when effort weight was 
decreased, the result was closer to the tracking data resulting in higher knee flexion, higher knee moment and 
higher activation in knee extensors during stance. However, the other variables were predicted worse. Weights 
in the objective function could be obtained from data using inverse optimal control instead of selecting them 
 empirically41,51. Third, it is not yet known which energy measure is minimized in human walking. Several studies 
proposed that metabolic energy is  minimized52–54, while others hypothesized that it is more likely that muscular 
effort, which is related to activation and thus neural excitation, is minimized in human  gait22,49,55. Since the actual 
movement objective is unknown, we included tracking. The simulation is still predictive because a new motion 
task was simulated based on another one. A different option is to manually tune weightings of energy measures 
to predict walking and  running24. However, instead of data, this requires expert input.
In conclusion, we presented a comprehensive 3D full-body musculoskeletal model modified for biomechani-
cal analysis of running with directional change, i.e. curved running. Model dynamics were formulated implicitly 
resulting in computational efficient simulations. The efficiency makes large scale inverse optimal control stud-
ies or sensitivity studies actually feasible. Furthermore, virtual product  design11,13 would considerably benefit. 
Predicted kinematics and kinetics confirmed the predictive power of the proposed approach and were very 
promising but limited by the fact that the true objective of human motion is still unknown. For this reason, this 
work might be an important step towards efficient and biomechanical accurate predictive simulations of move-
ments including directional changes.
Received: 30 April 2020; Accepted: 21 September 2020
References
 1. Ezati, M., Ghannadi, B. & McPhee, J. A review of simulation methods for human movement dynamics with emphasis on gait. 
Multibody Syst. Dyn. 47, 265–292 (2019).
 2. Dorschky, E., Nitschke, M., Seifer, A.-K., Van den Bogert, A. J. & Eskofier, B. M. Estimation of gait kinematics and kinetics from 
inertial sensor data using optimal control of musculoskeletal models. J. Biomech. 95, 109278 (2019).
 3. Erdemir, A., McLean, S., Herzog, W. & Van den Bogert, A. J. Model-based estimation of muscle forces exerted during movements. 
Clin. Biomech. 22, 131–154 (2007).
 4. Pandy, M. G. Computer modeling and simulation of human movement. Annu. Rev. Biomed. Eng. 3, 245–273 (2001).
 5. Seth, A. et al. Opensim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. 
PLoS Comput. Biol. 14, 1–20 (2018).
 6. Valero-Cuevas, F. J., Hoffmann, H., Kurse, M. U., Kutch, J. J. & Theodorou, E. A. Computational models for neuromuscular func-
tion. IEEE Rev. Biomed. Eng. 2, 110–135 (2009).
 7. Van den Bogert, A. J., Blana, D. & Heinrich, D. Implicit methods for efficient musculoskeletal simulation and optimal control. 
Procedia IUTAM 2, 297–316 (2011).
 8. Arnold, A. S., Liu, M. Q., Schwartz, M. H., Ounpuu, S. & Delp, S. The role of estimating muscle-tendon lengths and velocities of 
the hamstrings in the evaluation and treatment of crouch gait. Gait Posture 23, 273–281 (2006).
 9. Goldberg, S. R., Õunpuu, S. & Delp, S. L. The importance of swing-phase initial conditions in stiff-knee gait. J. Biomech. 36, 
1111–1116 (2003).
 10. Fey, N. P., Klute, G. K. & Neptune, R. R. Optimization of prosthetic foot stiffness to reduce metabolic cost and intact knee loading 
during below-knee amputee walking: A theoretical study. J. Biomech. Eng. 134, 1–10 (2012).
 11. Koelewijn, A. D. & Van den Bogert, A. J. Joint contact forces can be reduced by improving joint moment symmetry in below-knee 
amputee gait simulations. Gait Posture 49, 219–225 (2016).
 12. Harant, M., Sreenivasa, M., Millard, M., Šarabon, N. & Mombaur, K. Parameter optimization for passive spinal exoskeletons based 
on experimental data and optimal control. In 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids), 
535–540 (2017).
 13. Dorschky, E., Krüger, D., Kurfess, N., Schlarb, H., Wartzack, S., Eskofier, B. M. & Van den Bogert, A. J. Optimal control simulation 
predicts effects of midsole materials on energy cost of running. Comput. Methods Biomech. Biomed. Eng. 22, 869–879 (2019).
 14. Fluit, R., Andersen, M., Kolk, S., Verdonschot, N. & Koopman, H. Prediction of ground reaction forces and moments during 
various activities of daily living. J. Biomech. 47, 2321–2329 (2014).
 15. Geyer, H. & Herr, H. A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and 
muscle activities. IEEE Trans. Neural Syst. Rehabil. Eng. 18, 263–273 (2010).
 16. Song, S. & Geyer, H. A neural circuitry that emphasizes spinal feedback generates diverse behaviours of human locomotion. J. 
Physiol. 593, 3493–3511 (2015).
 17. Kim, Y., Tagawa, Y., Obinata, G. & Hase, K. Robust control of CPG-based 3D neuromusculoskeletal walking model. Biol. Cybern. 
105, 269–282 (2011).
 18. McLean, S. G., Huang, X. & Van den Bogert, A. J. Investigating isolated neuromuscular control contributions to non-contact 
anterior cruciate ligament injury risk via computer simulation methods. Clin. Biomech. 23, 926–936 (2008).
 19. Taboga, P. & Kram, R. Modelling the effect of curves on distance running performance. PeerJ 7, e8222 (2019).
 20. Song, S. & Geyer, H. Regulating speed and generating large speed transitions in a neuromuscular human walking model. In 2012 
IEEE International Conference on Robotics and Automation, 511–516 (2012).
 21. Miller, R. H. & Hamill, J. Computer simulation of the effects of shoe cushioning on internal and external loading during running 
impacts. Comput. Methods Biomech. Biomed. Eng. 12, 481–490 (2009).
 22. Miller, R. H., Umberger, B. R., Hamill, J. & Caldwell, G. E. Evaluation of the minimum energy hypothesis and other potential 
optimality criteria for human running. Proc. R. Soc. B Biol. Sci. 279, 1498–1505 (2012).
11
Vol.:(0123456789)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
 23. Miller, R. H. & Hamill, J. Optimal footfall patterns for cost minimization in running. J. Biomech. 48, 2858–2864 (2015).
 24. Falisse, A., Serrancoli, G., Dembia, C. L., Gillis, J., Jonkers, I., De Groote, F. Rapid predictive simulations with complex muscu-
loskeletal models suggest that diverse healthy and pathological human gaits can emerge from similar control strategies. J. R. Soc. 
Interface 16, 20190402 (2019).
 25. Lin, Y.-C. & Pandy, M. G. Three-dimensional data-tracking dynamic optimization simulations of human locomotion generated 
by direct collocation. J. Biomech. 59, 1–8 (2017).
 26. Lin, Y.-C., Walter, J. P. & Pandy, M. G. Predictive simulations of neuromuscular coordination and joint-contact loading in human 
gait. Ann. Biomed. Eng. 46(8), 1216–1227 (2018).
 27. Miller, R. H. A comparison of muscle energy models for simulating human walking in three dimensions. J. Biomech. 47, 1373–1381 
(2014).
 28. Lee, L.-F. & Umberger, B. R. Generating optimal control simulations of musculoskeletal movement using opensim and matlab. 
PeerJ 4, e1638 (2016).
 29. Mordatch, I., Mishra, N., Eppner, C. & Abbeel, P. Combining model-based policy search with online model learning for control of 
physical humanoids. In 2016 IEEE International Conference on Robotics and Automation (ICRA), vol. 2016-June, 242–248 (IEEE, 
2016).
 30. Porsa, S., Lin, Y.-C. & Pandy, M. G. Direct methods for predicting movement biomechanics based upon optimal control theory 
with implementation in opensim. Ann. Biomed. Eng. 44, 2542–2557 (2016).
 31. Delp, S. L. et al. Opensim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. 
Eng. 54, 1940–1950 (2007).
 32. Hamner, S., Seth, A. & Delp, S. L. Muscle contributions to propulsion and support during running. J. Biomech. 43, 2709–2716 
(2010).
 33. Baker, R. Pelvic angles: a mathematically rigorous definition which is consistent with a conventional clinical understanding of the 
terms. Gait Posture 13, 1–6 (2001).
 34. Hill, A. V. The heat of shortening and the dynamic constants of muscle. Proc. R. Soc. Lond. B Biol. Sci. 126, 136–195 (1938).
 35. Kyröläinen, H., Komi, P. V. & Belli, A. Changes in muscle activity patterns and kinetics with increasing running speed. J. Strength 
Cond. Res. 13, 400 (1999).
 36. Neptune, R. R., Sasaki, K. & Kautz, S. A. The effect of walking speed on muscle function and mechanical energetics. Gait Posture 
28, 135–143 (2008).
 37. Van den Bogert, A. J. & de Koning, J. On optimal filtering for inverse dynamics analysis. In Proceedings of the IXth Biennial Confer-
ence of the Canadian Society for Biomechanics., 214–215 (Vancouver, British Columbia, 1996).
 38. Derrick, T. R. et al. ISB recommendations on the reporting of intersegmental forces and moments during human motion analysis. 
J. Biomech. 99, 109533 (2020).
 39. Wächter, A. & Biegler, L. T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear 
programming. Math. Program. 106, 25–57 (2006).
 40. Cappellini, G., Ivanenko, Y. P., Poppele, R. E. & Lacquaniti, F. Motor patterns in human walking and running. J. Neurophysiol. 95, 
3426–3437 (2006).
 41. Clever, D., Malin Schemschat, R., Felis, M. L. & Mombaur, K. Inverse optimal control based identification of optimality criteria in 
whole-body human walking on level ground. In 2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics 
(BioRob), 1192–1199 (IEEE, 2016).
 42. Felis, M. L. & Mombaur, K. Synthesis of full-body 3-D human gait using optimal control methods. In 2016 IEEE International 
Conference on Robotics and Automation (ICRA), vol. 2016-June, 1560–1566 (IEEE, 2016).
 43. Falisse, A., Serrancolí, G., Dembia, C. L., Gillis, J. & De Groote, F. Algorithmic differentiation improves the computational efficiency 
of opensim-based trajectory optimization of human movement. PLoS ONE 14, e0217730 (2019).
 44. Jiang, Y., Van Wouwe, T., De Groote, F. & Liu, C. K. Synthesis of biologically realistic human motion using joint torque actuation. 
ACM Trans. Graph. 38, 1–12 (2019).
 45. Lee, S., Lee, K., Park, M. & Lee, J. Scalable muscle-actuated human simulation and control. ACM Trans. Graph. 38(4), 1–13 (2019).
 46. Akrami, M. et al. Subject-specific finite element modelling of the human foot complex during walking: sensitivity analysis of 
material properties, boundary and loading conditions. Biomech. Model. Mechanobiol. 17, 559–576 (2018).
 47. Bogo, F. et al. Keep it SMPL: Automatic estimation of 3D human pose and shape from a single image. In Computer Vision—ECCV 
2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part V, 561–578 (Springer Inter-
national Publishing, 2016).
 48. Koelewijn, A. D. & Van den Bogert, A. J. A solution method for predictive simulations in a stochastic environment. J. Biomech. 
104, 109759. https ://doi.org/10.1016/j.jbiom ech.2020.10975 9 (2020).
 49. Ackermann, M. & Van den Bogert, A. J. Optimality principles for model-based prediction of human gait. J. Biomech. 43, 1055–1060 
(2010).
 50. De Groote, F., Kinney, A. L., Rao, A. V. & Fregly, B. J. Evaluation of direct collocation optimal control problem formulations for 
solving the muscle redundancy problem. Ann. Biomed. Eng. 44, 2922–2936 (2016).
 51. Nguyen, V. Q., Johnson, R. T., Sup, F. C. & Umberger, B. R. Bilevel optimization for cost function determination in dynamic 
simulation of human gait. IEEE Trans. Neural Syst. Rehabil. Eng. 27, 1426–1435 (2019).
 52. Gordon, K. E., Ferris, D. P. & Kuo, A. D. Metabolic and mechanical energy costs of reducing vertical center of mass movement 
during gait. Arch. Phys. Med. Rehabil. 90, 136–144 (2009).
 53. Ralston, H. J. Energy-speed relation and optimal speed during level walking. Internationale Zeitschrift für Angewandte Physiologie 
Einschliesslich Arbeitsphysiologie 17, 277–283 (1958).
 54. Zarrugh, M. Y., Todd, F. N. & Ralston, H. J. Optimization of energy expenditure during level walking. Eur. J. Appl. Physiol. 33, 
293–306 (1974).
 55. McDonald, K. A., Hieronymi, A., Cusumano, J. P. & Rubenson, J. Optimization in human walking: Decoupling whole-body ener-
getics and local muscle effort. In ISB/ASB Conference 2019 (2019).
Acknowledgements
This work was supported by adidas AG (M.N., A.D.K., A.J.v.d.B.), the Bavarian Ministry of Economic Affairs, 
Regional Development and Energy within the Embedded Systems Initiative (E.D.), the German Research Foun-
dation within the framework of the Heisenberg professorship programme under Grant ES 434/8-1 (B.M.E.), the 
national science foundation under Grant No. 1344954 (A.D.K., A.J.v.d.B.), and a Graduate Scholarship from the 
Parker-Hannifin Corporation (A.D.K.).
Author contributions
A.J.v.d.B. had the original idea of using an implicit formulation of the model dynamics. M.N., E.D., A.D.K, and 
A.J.v.d.B. developed and implemented the model, which was supported by D.H., H.S., and B.M.E. The data 
12
Vol:.(1234567890)
Scientific Reports |        (2020) 10:17655  | 
https://doi.org/10.1038/s41598-020-73856-w
www.nature.com/scientificreports/
was recorded by H.S. The simulations and the manuscript were created by M.N. All authors contributed to the 
interpretation of the results and to the editing of the manuscript.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Competing interests 
The authors declare no competing interests.
Additional information
Supplementary information is available for this paper at https ://doi.org/10.1038/s4159 8-020-73856 -w.
Correspondence and requests for materials should be addressed to M.N.
Reprints and permissions information is available at www.nature.com/reprints.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and 
institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International 
License, which permits use, sharing, adaptation, distribution and reproduction in any medium or 
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the 
Creative Commons licence, and indicate if changes were made. The images or other third party material in this 
article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the 
material. If material is not included in the article’s Creative Commons licence and your intended use is not 
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.
© The Author(s) 2020

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