# Change the direction: 3D optimal control simulation by directly tracking marker and ground reaction force data

- Published
- Accepted
- Received

- Academic Editor
- Mike Climstein

- Subject Areas
- Computational Biology, Kinesiology, Radiology and Medical Imaging, Computational Science, Biomechanics
- Keywords
- Dynamic optimization, Trajectory optimization, Human Movement Simulation, Motion capturing, Musculoskeletal model, Locomotion, Three-dimensional

- Copyright
- © 2023 Nitschke et al.
- Licence
- This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.

- Cite this article
- 2023. Change the direction: 3D optimal control simulation by directly tracking marker and ground reaction force data. PeerJ 11:e14852 https://doi.org/10.7717/peerj.14852

## Abstract

Optimal control simulations of musculoskeletal models can be used to reconstruct motions measured with optical motion capture to estimate joint and muscle kinematics and kinetics. These simulations are mutually and dynamically consistent, in contrast to traditional inverse methods. Commonly, optimal control simulations are generated by tracking generalized coordinates in combination with ground reaction forces. The generalized coordinates are estimated from marker positions using, for example, inverse kinematics. Hence, inaccuracies in the estimated coordinates are tracked in the simulation. We developed an approach to reconstruct arbitrary motions, such as change of direction motions, using optimal control simulations of 3D full-body musculoskeletal models by directly tracking marker and ground reaction force data. For evaluation, we recorded three trials each of straight running, curved running, and a v-cut for 10 participants. We reconstructed the recordings with marker tracking simulations, coordinate tracking simulations, and inverse kinematics and dynamics. First, we analyzed the convergence of the simulations and found that the wall time increased three to four times when using marker tracking compared to coordinate tracking. Then, we compared the marker trajectories, ground reaction forces, pelvis translations, joint angles, and joint moments between the three reconstruction methods. Root mean squared deviations between measured and estimated marker positions were smallest for inverse kinematics (*e.g*., 7.6 ± 5.1 mm for v-cut). However, measurement noise and soft tissue artifacts are likely also tracked in inverse kinematics, meaning that this approach does not reflect a gold standard. Marker tracking simulations resulted in slightly higher root mean squared marker deviations (*e.g*., 9.5 ± 6.2 mm for v-cut) than inverse kinematics. In contrast, coordinate tracking resulted in deviations that were nearly twice as high (*e.g*., 16.8 ± 10.5 mm for v-cut). Joint angles from coordinate tracking followed the estimated joint angles from inverse kinematics more closely than marker tracking (*e.g*., root mean squared deviation of 1.4 ± 1.8 deg *vs*. 3.5 ± 4.0 deg for v-cut). However, we did not have a gold standard measurement of the joint angles, so it is unknown if this larger deviation means the solution is less accurate. In conclusion, we showed that optimal control simulations of change of direction running motions can be created by tracking marker and ground reaction force data. Marker tracking considerably improved marker accuracy compared to coordinate tracking. Therefore, we recommend reconstructing movements by directly tracking marker data in the optimal control simulation when precise marker tracking is required.

## Introduction

Kinematics and kinetics of walking or running are estimated from measurements with optical motion capture in various fields of biomechanical research. While research often focuses on straight walking or running, change of direction (COD) motions are also crucial in everyday life. Cutting maneuvers are, for example, performed frequently in multi-directional team sports (Fox, 2018). Non-contact COD maneuvers or rapid decelerations have been identified as the primary cause of anterior cruciate ligament (ACL) injuries (Donnelly et al., 2017; McLean et al., 2004). In the last years, a great number of biomechanical studies (*e.g*., Barengo et al. (2014)) were conducted to develop and analyze the effect of injury prevention training programs like FIFA 11+ (Bizzini, Junge & Dvorak, 2011). However, in those studies, little attention was spent on the method that was used to estimate the kinematic and kinetic variables like joint angles, joint moments, or muscle forces.

Inverse methods, *i.e*., inverse kinematics and dynamics, combined with static or dynamic optimization using a human model are widely used to obtain joint and muscle kinematics and kinetics from marker positions and ground reaction forces (GRFs) (Seth et al. (2018); see Fig. 1A). However, whereas these methods have the advantage that they are rapid to solve and easy to apply, they also have major weaknesses. Inverse kinematics estimates the generalized coordinates of the model, *i.e*., global translation, global orientation, and joint angles, for each time step separately. Since time dependency is not taken into account, inverse kinematics is prone to track measurement noise. This means that, when measurement noise causes a sudden change in a marker position, inverse kinematics will also estimate a sudden change in the respective joint angle, even though it is unrealistic for humans to move in a non-smooth fashion. In a second step, joint moments are estimated with inverse dynamics. While inverse dynamics takes time dependency into account, it allows for dynamic inconsistencies, *i.e*., inconsistencies between kinematics and kinetics. These inconsistencies are caused by modeling errors and inaccuracies in the measured data and typically result in residual forces and moments at the last segment (Faber, van Soest & Kistemaker, 2018). However, there is no physical cause for these residuals and it is difficult to trace which model parameter or measurement error contributed to the residuals. Therefore, residuals are hard to interpret and researchers are advised to prevent residuals that are large enough to influence the study conclusion (Hicks et al., 2015). After inverse dynamics, muscle forces can be computed using static or dynamic optimization to resolve the muscle redundancy problem. Past studies found minor differences in muscle forces obtained from dynamic compared to static optimization for walking (Anderson & Pandy, 2001) and even for running (Lin et al., 2012). However, more recent evidence highlights that muscle activation, and thus muscle efficiency, is influenced when neglecting tendon compliance (De Groote et al., 2016; Miller, Umberger & Caldwell, 2012). Hence, modeling muscle dynamics is especially important for faster motions such as running or sprinting.

Open-loop optimal control simulation of a human model is an alternative to inverse methods for estimating various biomechanical variables from a measured motion, which results in mutually and dynamically consistent kinematics and kinetics. In open-loop optimal control simulations, also called trajectory optimization, joint and muscle kinematics and kinetics of the model are obtained by minimizing an objective while accounting for system dynamics. When open-loop optimal control simulations are used to reconstruct movements, it is assumed that the skeletal or musculoskeletal model is perfect and does not contain modeling errors. For reconstructing a measured movement with optimal control simulations, the objective often combines a tracking term minimizing the difference between measured and simulated data with an energy-related term. Optimal control simulations have gained increased attention in recent years due to methodological advances. The exploration of direct collocation methods combined with the implicit formulation of the system dynamics allows the optimal control problem to be solved efficiently (De Groote et al., 2016; Nitschke et al., 2020; van den Bogert, Blana & Heinrich, 2011). Furthermore, the increasing availability of toolboxes facilitates access to the methodology (Dembia et al., 2020; Michaud et al., 2022; Patterson & Rao, 2014). Optimal control simulations cannot only be used to reconstruct measured motions but also to predict responses to environmental changes (Dorschky et al., 2019; van den Bogert, Blana & Heinrich, 2011; van den Bogert et al., 2012) or task changes (Lin, Walter & Pandy, 2018; Nitschke et al., 2020). Predictive simulations can either track related data as reference or predict novel movements without any input from measurements.

Reconstructive optimal control simulations are traditionally generated by tracking generalized coordinates or joint angles in combination with measured GRFs (*e.g*., Dembia et al. (2020); Haralabidis et al. (2021); Heinrich, van den Bogert & Nachbauer (2014); Lin, Walter & Pandy (2018); Nitschke et al. (2020); van den Bogert, Blana & Heinrich (2011); van den Bogert et al. (2012); see Fig. 1B). Since generalized coordinates are used as kinematic states of the model and are thus optimization variables, tracking of coordinates is computationally more efficient than tracking of other biomechanical variables which are not part of the optimization variables. However, the coordinates have to be estimated from marker data before simulation using, for example, inverse kinematics. Hence, the inaccuracies in the estimated coordinates are tracked in the simulation resulting in error propagation. Additionally, individual joint angles are tracked rather than absolute positions. Therefore, tracking errors of each joint angle accumulate down the kinematic chain, which can cause larger positional differences at the end of this chain. In contrast to tracking coordinates in the simulation, tracking marker positions directly could avoid error propagation and error accumulation along the kinematic chain (see Fig. 1C). Recently, marker tracking was successfully investigated for upper limb models with up to 7 degrees of freedom (DoFs) and up to 20 muscle tendon units (MTUs) partly in combination with electromyography (EMG) tracking (Bailly et al., 2021; Bélaise et al., 2018a, 2018b; Hoffmann et al., 2020). Furthermore, marker tracking was investigated for a single leg model with 6 DoFs and 17 MTUs (Moissenet et al., 2019). Febrer-Nafría et al. (2020) and Venne et al. (2022) compared coordinate and marker tracking using full-body skeletal models for walking and somersaults, respectively. Their research indicated that marker tracking simulations followed measured marker positions more closely than coordinate tracking. Therefore, marker tracking was more accurate in terms of marker errors.

Overall, previous research on marker tracking was limited to small models with few optimization variables and was limited to an evaluation with simulated or little data. The models previously used for marker tracking were either musculoskeletal models with only a few DoFs (Bailly et al., 2021; Bélaise et al., 2018a, 2018b; Moissenet et al., 2019) or skeletal models (Febrer-Nafría et al., 2020; Hoffmann et al., 2020; Venne et al., 2022). Therefore, it is unclear whether optimal control simulation with marker tracking is numerically feasible for full-body musculoskeletal models, resulting in a considerably larger number of optimization variables and constraints. However, entire 3D body kinematics and kinetics should be considered especially for an accurate analysis of COD running motions since upper-body kinematics can influence, for example, knee moments (Donnelly et al., 2012). Moreover, COD running motions have not yet been reconstructed with optimal control simulation, which would particularly be relevant for sports science. Furthermore, evaluation was performed either with simulated data (Bailly et al., 2021; Bélaise et al., 2018a) or data of only one participant (Bélaise et al., 2018b; Febrer-Nafría et al., 2020; Hoffmann et al., 2020; Moissenet et al., 2019) except for Venne et al. (2022), who reconstructed in total 26 somersault trials of five participants. Consequently, there is no clear evidence of whether marker tracking is superior to coordinate tracking for running and especially for COD running motions.

We investigated the feasibility of directly driving 3D optimal control simulations by marker and GRF data using a full-body musculoskeletal model, especially for reconstructing COD running motions. We developed a method for creating optimal control simulations for arbitrary motions like COD motions with direct collocation and an implicit formulation of the system dynamics. Since gold standard measurements of kinematics are hardly available and joint kinetics cannot be measured directly, we compared marker tracking simulations to coordinate tracking simulations and inverse methods for estimated marker positions, GRFs, pelvis translation, angles, and joint moments (see Fig. 1). To create strong evidence with our study, we performed the analysis for 10 participants and three trials each of straight running, curved running, and a v-cut.

## Methods

In this section, we first describe the experimental data and the musculoskeletal model used for motion reconstruction. Then, we give details about the inverse methods, the optimal control simulations, and the evaluation.

### Experimental data

We recorded motion capture data of 10 healthy young participants (four female, six male; age: 27.5 $\pm $ 3.5 years; height: 1.76 $\pm $ 0.10 m; mass: 71.3 $\pm $ 12.1 kg). The ethics committee of the Friedrich-Alexander-Universität Erlangen-Nürnberg (Re.-No. 106_13 B) approved the study, and participants gave informed written consent before participation. We obtained marker positions of 42 reflective markers at 175 Hz with 11 infrared cameras (Qualisys, Gothenburg, Sweden) and GRFs of the right and left foot at 1750 Hz with two force plates (Bertec Corporation, Columbus, USA). Simultaneously, we recorded data from 11 inertial measurement units but did not use the data in this article.

Each participant first performed a static trial in a neutral pose (N-pose) with one foot on each force plate and the arms beside the body. Afterwards, the participants completed multiple trials respectively for straight running, curved running with a radius of 7 m, and a 90° v-cut (see Fig. 2). For curved running and the v-cut, we indicated the path with crepe tape on the floor.

For every participant and motion type, we choose the three trials for which the force plates were hit entirely and marker occlusions were smallest. We filled gaps in the marker data using the Qualisys Track Manager but did not apply any filter. We defined, but did not extract, the motion of interest before reconstruction to reduce edge effects from filtering for inverse dynamics or from the initial constraint used in the simulation. The motions of interest started at the initial contact of the right foot at the first force plate and ended with the next initial contact of the right foot. For the v-cut, this corresponded to the execution and departure contact. We determined each initial contact at the time where the vertical velocity of the mean of the heel and toe marker position was at a minimum (O’Connor et al., 2007).

### Musculoskeletal model

We used our 3D full-body musculoskeletal model called runMaD, which is short for “running model for motions in all directions” (Nitschke et al., 2020). In contrast to other musculoskeletal models, runMaD has an adapted pelvis rotation sequence that makes the pelvis obliquity and tilt interpretable according to their clinical definition independent of the movement direction, *i.e*., independent of the rotation around the vertical axis (Baker, 2001). The model has 33 DoFs (six DoFs between ground and pelvis, three DoFs at the lumbar joint, seven DoFs per leg, five DoFs per arm), 92 MTUs in the lower body (six MTUs at the lumbar joint,
$43$ MTUs per leg), and five torque actuators per arm. The muscle paths are described in OpenSim using point sets. For the optimal control simulation, muscle-tendon lengths are defined with polynomial functions depending on the joint angles, and a penetration-based ground contact model with eight contact points at each foot is used (see supplemental information of Nitschke et al., 2020). We scaled the model for every participant using the static trial in N-pose in OpenSim 4.3 (Seth et al., 2018) but did not personalize any muscle properties.

### Inverse methods

Using the scaled models, we performed inverse kinematics and dynamics as reference with OpenSim 4.3 (Seth et al., 2018). Additionally, we used the resulting generalized coordinates as input for the coordinate tracking simulations (see Fig. 1). We weighted all markers equally in inverse kinematics. For inverse dynamics, we filtered the generalized coordinates and GRFs with a 3rd order dual-pass low-pass Butterworth filter with a cut-off frequency of 15 Hz (Derrick et al., 2020). We reconstructed the entire trials and not only the extracted motions of interest since the Butterworth filter has an infinite impulse response causing undesired effects, especially at the edges of the trajectories.

### Optimal control simulations

We created coordinate and marker tracking simulations by solving optimal control problems with the scaled musculoskeletal models. We formulated the optimal control problem as a constrained non-linear optimization problem using direct collocation and a backward Euler discretization. A state trajectory $\mathbf{\text{x}}$ and a control trajectory $\mathbf{\text{u}}$ of the model are found by minimizing a multi-objective function $J(\mathbf{x}\mathbf{,}\mathbf{u}\mathbf{)}$ with respect to the model dynamics $\mathbf{\text{f}}$.

**Objective Function** The objective *J* was a weighted sum of tracking
${J}_{tra}$, muscular effort
${J}_{mus}$, torque effort
${J}_{tor}$, and regularization
${J}_{reg}$:

(1) $$J(\mathbf{x}\mathbf{,}\mathbf{u}\mathbf{)}\mathbf{=}{\mathbf{J}}_{\mathbf{t}\mathbf{r}\mathbf{a}}\mathbf{+}{\mathbf{J}}_{\mathbf{m}\mathbf{u}\mathbf{s}}\mathbf{+}{\mathbf{J}}_{\mathbf{t}\mathbf{o}\mathbf{r}}\mathbf{+}{\mathbf{J}}_{\mathbf{r}\mathbf{e}\mathbf{g}}\mathbf{,}$$

(2) $${J}_{tra}=\sum _{j\in {S}_{tra}}\frac{{W}_{tra,j}}{N{N}_{tra,j}}\sum _{k=1}^{N}\sum _{i=1}^{{N}_{tra,j}}{\left({y}_{i,j}\right[k]-{y^}_{i,j}[k\left]\right)}^{2}$$

(3) $${J}_{mus}=\frac{{W}_{mus}}{N{N}_{mus}}\text{}\sum _{k=1}^{N}\text{}\sum _{i=1}^{{N}_{mus}}\frac{{w}_{mus,i}}{{\displaystyle \sum _{j=1}^{{N}_{mus}}}{w}_{mus,j}}{\left({n}_{e,i}\right[k\left]\right)}^{3},$$

(4) $${J}_{tor}=\frac{{W}_{tor}}{N{N}_{tor}}\sum _{k=1}^{N}\sum _{i=1}^{{N}_{tor}}{\left({m}_{i}\right[k\left]\right)}^{2},$$

(5) $$\begin{array}{ll}& {J}_{reg}=\frac{{W}_{reg}(N-1)}{{T}^{2}({N}_{states}+{N}_{controls})}\text{}\sum _{k=1}^{N-1}\\ & {\textstyle \phantom{\rule{2em}{0ex}}}(\sum _{i=1}^{{N}_{states}}({x}_{i}[k+1]-{x}_{i}[k]{)}^{2}+\sum _{i=1}^{{N}_{controls}}({u}_{i}[k+1]-{u}_{i}[k]{)}^{2}).\end{array}$$The tracking term ${J}_{tra}$ consisted of separate terms with individual weights ${W}_{tra,j}$ for each data type $j$ of the set ${S}_{tra}$. Depending on the tracking method, we used the following data types (see Fig. 1):

**Coordinate tracking:**3D global translation of the pelvis, global orientation of the pelvis, and joint angles obtained from inverse kinematics and measured GRFs of right and left foot (*i.e*., ${S}_{tra}=\{translation,angle,\text{GRF}\}$, ${N}_{tra,translation}=3$, ${N}_{tra,angle}=30$, and ${N}_{tra,GRF}=6$).**Marker tracking:**measured 3D marker positions of all 42 markers and GRFs of right and left foot (*i.e*., ${S}_{tra}=\{marker,\text{GRF}\}$, ${N}_{tra,marker}=126$, and ${N}_{tra,GRF}=6$).

We minimized the squared difference between reference signal
$y$ and estimated signal
$\hat{y}$ in the tracking term for *N* collocation nodes and
${N}_{tra,j}$ signals. Using
${J}_{mus}$ with the respective weight
${W}_{mus}$, we resolved the muscle redundancy problem by minimizing the sum of the volume-weighted cubed neural excitations
${n}_{e}$ of each of the
${N}_{mus}$ muscles. We used the muscle volume
${w}_{mus,i}$ of a muscle
$i$ to account in the effort term for the strongly varying sizes and maximum isometric forces of the MTUs and, therefore, spread muscle recruitment more evenly (Happee & Van der Helm, 1995).

Furthermore, we minimized the sum of the squared torque controls
$m$ actuating the
${N}_{tor}$ arms in
${J}_{tor}$ with the weight
${W}_{tor}$. The regularization term
${J}_{reg}$ with the weight
${W}_{reg}$ represents the minimization of the temporal derivative of the state and control trajectories, where
${N}_{states}$ and
${N}_{controls}$ represent the number of states and controls, respectively. Such regularizations are used to improve the convergence of the optimization algorithm. The duration *T* of the motion was prescribed by the duration of the tracking data.

**Model Dynamics** We formulated the model dynamics implicitly and used backward Euler discretization. Hence, the following constraint was applied for each node:

(6) $$\mathbf{f}(\mathbf{x}[k+1],\frac{\mathbf{x}[k+1]-\mathbf{x}[k]}{h},\mathbf{u}[k+1])=\mathbf{0}\phantom{\rule{1em}{0ex}}\mathrm{\forall}k=1,...,N-1\phantom{\rule{thinmathspace}{0ex}},$$where $h=T/(N-1)$. More details about the system dynamics and the implementation are given by Nitschke et al. (2020).

However, Eq. (6) does not contain the control $\mathbf{u}[1]$, as the following dynamics apply at the first node $k=1$:

(7)
$$\mathbf{f}(\mathbf{x}[2],\frac{\mathbf{x}[2]-\mathbf{x}[1]}{h},\mathbf{u}[2])=\mathbf{0}.$$Furthermore, the state
$\mathbf{x}[1]$ at the first node appears only in one equation (equation
$k=1$), while all other states
$\mathbf{x}[k]$ appear in two subsequent equations (equation
$k-1$ and equation
$k$). Hence, additional information is required to ensure that the optimal control problem can be solved, which is commonly done in two ways. The first option is to add an initial state. The second option is to apply an additional constraint that describes the task. For example, gait is typically constrained to be periodic (*e.g*., van den Bogert, Blana & Heinrich (2011)). The usage of task constraints can especially be beneficial when a new motion should be predicted based on data of a related motion (Nitschke et al., 2020) or when a specific gait speed should be prescribed for standardization (Dorschky et al., 2019). In this work, however, we aimed to reconstruct arbitrary motions, which would not be possible if an initial state or task constraint had to be prescribed. Therefore, we did not use a task constraint but additionally ensured model dynamics using forward Euler discretization at
$k=1$:

(8) $$\mathbf{f}(\mathbf{x}[1],\frac{\mathbf{x}[2]-\mathbf{x}[1]}{h},\mathbf{u}[1])=\mathbf{0}\phantom{\rule{thinmathspace}{0ex}},$$except for the identities $\dot{\mathbf{q}}-\frac{d\mathbf{q}}{\frac{}{}}dt=\mathbf{0}$ of the global pelvis translation and orientation to not restrict global motion. The combination of Eqs. (6) and (8) implies constant velocities of the states and controls between nodes 1 and 2. In contrast to prescribing specific values as initial states, this constraint does not require prior knowledge of the motion. However, the assumption of a constant velocity slightly influences the result at the first nodes. To avoid impact on the motion of interest, we included additional samples at the beginning of the signal.

**Initialization** We first simulated static standing for each participant and each tracking method to calibrate the ground contact model and to generate an initial guess for the running simulations. In the objective, we tracked the kinematic and GRF data of one time point of the N-pose. Since only one time point was simulated, we omitted the regularization term
${J}_{reg}$ and ensured static equilibrium by
$\mathbf{f}(\mathbf{x}[1],\mathbf{0},\mathbf{u}[1])=\mathbf{0}$, which ensures that the velocities and accelerations are zero. We adapted the ground contact model for different shoe sole thicknesses of the participants by optimizing a vertical offset for the position of the ground contact points during the simulation. The position of the ground contact points had to be adapted for coordinate and marker tracking since we tracked absolute positions of the pelvis or markers, respectively. We solved 10 simulations for each optimization problem using different random initial guesses and selected the solution with the lowest objective to reduce the chance of obtaining a local minimum.

**Running Simulations** We then generated simulations for three trials each of straight running, curved running, and the v-cut for each of the 10 participants and tracking method. In total, this resulted in 90 simulations each for marker and coordinate tracking. In the objective, we tracked reference kinematics and GRFs using the sampling frequency of 175 Hz and therefore downsampled the GRF data. The reference data was not filtered prior to the simulation since the simulation itself acts as a physical filter that takes the model dynamics into account. We found in pilot simulations that 10 additional samples before the motion of interest, *i.e*., before the initial contact, are sufficient to not cause observable artifacts in the motion of interest when using Eq. (8). The original sampling frequency combined with the additional samples resulted in optimization problems with 127 to 180 collocation nodes *N*, which corresponds to durations *T* of approximately 0.72 s to 1.02 s.

**Solution Process** We selected the weights *W* of the objective terms (see Eqs. (2)–(5)) empirically using data of the first participant such that the tracking data was followed while the neural excitation remained smooth (see Table 1). We used an equal weight for all tracking variables of the same type, meaning that we used one weight for pelvis translations, one for angles, one for markers, and one for GRFs. The constrained non-linear optimization problems were solved using IPOPT 3.12.3 (Wächter & Biegler, 2006) with a convergence tolerance for the scaled nonlinear program (NLP) error of
${10}^{-4}$ and a maximum number of iterations of
$2\cdot {10}^{4}$. We used a high-performance cluster to parallelize the 180 running simulations. Each simulation was performed on a single cluster node with one Xeon E3-1240 CPU with four cores.

Standing | Running | |||
---|---|---|---|---|

Coordinate | Marker | Coordinate | Marker | |

${W}_{tra,translation}$ in mm ${}^{-2}$ | ${10}^{-3}$ | – | ${10}^{-3}$ | – |

${W}_{tra,angle}$ in deg ${}^{-2}$ | ${10}^{-1}$ | – | ${10}^{-1}$ | – |

${W}_{tra,marker}$ in mm ${}^{-2}$ | – | ${10}^{-2}$ | – | ${10}^{-2}$ |

${W}_{tra,GRF}$ in (BW%) ${}^{-2}$ | ${10}^{-2}$ | ${10}^{-2}$ | ${10}^{-3}$ | ${10}^{-3}$ |

${W}_{mus}$ | $1$ | $1$ | $1$ | $1$ |

${W}_{tor}$ | ${10}^{-1}$ | ${10}^{-1}$ | ${10}^{-1}$ | ${10}^{-1}$ |

${W}_{reg}$ | – | – | ${10}^{-3}$ | ${10}^{-3}$ |

### Evaluation

We analyzed the convergence of coordinate and marker tracking problems by comparing the number of iterations and the wall and CPU time required to solve the running simulations. Whereas the wall time represents the time passed to solve the problem, the CPU time captures the total time the single cores are active, *i.e.*, the total time required for calculations. The CPU time can therefore be greater than the wall time if multiple cores of a CPU are used for processing. Furthermore, we computed the CPU time per iteration to analyze the computational demand of a single iteration. To evaluate the computational demand of the objective, constraints, and their derivatives compared to the time spent in the optimization algorithm of IPOPT, we obtained the ratio of CPU time spent in the function evaluations of the NLP from the log file of IPOPT.

To investigate the estimated kinematics and kinetics, we extracted the motions of interest from the reconstructed motions. We visually compared the reconstructed motions of interest of inverse kinematics, coordinate tracking, and marker tracking using the visualization of the kinematics in OpenSim and trajectory graphs of marker positions, GRFs, pelvis translation, angles, and joint moments. Additionally, we obtained the root mean squared deviation (RMSD) for each motion of interest to analyze the agreement of the trajectories. We used the measured data as reference for marker positions and GRFs. For the generalized coordinates and joint moments, we considered the result of the inverse methods as reference since no ground truth was available. GRFs were scaled to body-weight percent (BW%) and joint moments were scaled to body-weight body-height percent (BW BH%). We aggregated all results by computing the mean and standard deviation over all variables of one type (*e.g*., marker positions) and all trials of one motion (*e.g*., straight running). Consequently, each mean value resulted from 30 simulations.

Finally, we evaluated residual forces and moments that are caused by dynamic inconsistencies. For inverse methods, we computed the root mean squared (RMS) residual forces and moments at the pelvis. In the coordinate and marker tracking simulations, we constrained the dynamic residuals to be zero by using the multibody dynamics as constraints in the optimization (see equation S1 in the supplemental information of Nitschke et al. (2020)). Nevertheless, the optimization result could slightly violate the multibody dynamics within the constraint violation tolerance of 0.001. Therefore, we analyzed the RMS residual pelvis forces, pelvis moments, and joint moments resulting from the constraint violations of the multibody dynamics. We scaled the residual forces to percent of maximal net ground reaction forces (GRF ${}_{\text{max}}$%) and the residual moments to percent of maximal net ground reaction forces and body-height percent (GRF ${}_{\text{max}}$ BH%) based on Hicks et al. (2015).

## Results

All 180 running simulations converged. The averages of the scaled NLP errors were between $7.3\cdot {10}^{-5}$ and $8.0\cdot {10}^{-5}$ for straight running, curved running, and v-cut and for marker and coordinate tracking (see Table 2). The NLP errors did not differ largely depending on the motion type or tracking method. Solving the marker tracking simulations required considerably more iterations and time than the coordinate tracking simulations. The fastest simulations for marker and coordinate tracking converged after a wall time of 2 h 42 min and 51 min, respectively, while the slowest simulations required 9 h 31 min and 3 h 56 min, respectively. The CPU time spent in every iteration was comparable between both tracking methods. However, the ratio of CPU time which was spent in the function evaluation of the NLP was higher for the marker tracking (10.1% to 10.5%) compared to the coordinate tracking (7.8% to 8.1%).

Motion | Marker tracking | Coordinate tracking | |
---|---|---|---|

Scaled NLP error | SR | 7.7 $\pm $ 1.5 $\cdot {10}^{-5}$ | 7.5 $\pm $ 1.7 $\cdot {10}^{-5}$ |

CR | 7.3 $\pm $ 2.0 $\cdot {10}^{-5}$ | 8.0 $\pm $ 1.6 $\cdot {10}^{-5}$ | |

VC | 7.3 $\pm $ 2.3 $\cdot {10}^{-5}$ | 7.7 $\pm $ 1.5 $\cdot {10}^{-5}$ | |

Number of iterations | SR | 7437 $\pm $ 2217 | 2204 $\pm $ 664 |

CR | 8068 $\pm $ 2030 | 2054 $\pm $ 505 | |

VC | 5788 $\pm $ 1699 | 2370 $\pm $ 821 | |

Wall time in hh:mm:ss | SR | 05:51:30 $\pm $ 01:48:54 | 01:45:42 $\pm $ 00:41:28 |

CR | 06:21:01 $\pm $ 01:27:17 | 01:34:10 $\pm $ 00:25:10 | |

VC | 04:57:36 $\pm $ 01:25:08 | 02:00:25 $\pm $ 00:39:40 | |

CPU time in hh:mm:ss | SR | 12:33:17 $\pm $ 03:47:35 | 03:54:00 $\pm $ 01:31:40 |

CR | 13:40:55 $\pm $ 03:02:47 | 03:30:17 $\pm $ 00:56:32 | |

VC | 10:40:31 $\pm $ 02:57:06 | 04:25:22 $\pm $ 01:27:08 | |

CPU time per iteration in ss:fff | SR | 06:115 $\pm $ 01:004 | 06:279 $\pm $ 00:763 |

CR | 06:196 $\pm $ 00:746 | 06:145 $\pm $ 00:802 | |

VC | 06:725 $\pm $ 01:027 | 06:794 $\pm $ 01:014 | |

CPU time in NLP in % | SR | 10.5 $\pm $ 1.6 | 7.8 $\pm $ 0.8 |

CR | 10.4 $\pm $ 1.2 | 8.1 $\pm $ 1.0 | |

VC | 10.1 $\pm $ 1.4 | 7.9 $\pm $ 1.1 |

The visual inspection of the reconstruction showed that the methods generally led to natural running motions. We overlaid the animated skeletons for inverse methods, coordinate tracking, and marker tracking for a more detailed kinematic analysis. Figure 3A shows an exemplary v-cut (participant 02, trial 130 in Nitschke et al. (2022)), also provided as a video in the supplemental information. The three reconstruction methods showed good agreement. Nevertheless, it could be observed that the result of the coordinate tracking deviated from that of inverse methods and marker tracking while the skeletons from inverse methods and marker tracking superimposed better. This can for example be seen for the right foot in the screenshots of the samples 31 and 121 in Fig. 3A.

For the different reconstruction methods, trajectories of marker positions, GRFs, joint angles, and joint moments showed similar patterns, but there were also substantial differences (see Fig. 3B). Mean RMSDs were in the same order of magnitude, except for marker positions (see Table 3). The mean RMSDs between estimated and measured marker positions were smallest for inverse methods (*e.g*., 7.6
$\pm $ 5.1 mm for v-cut) and a bit higher for the marker tracking simulation (*e.g*., 9.5
$\pm $ 6.2 mm for v-cut). However, coordinate tracking resulted in nearly twice as high RMSDs (*e.g*., 16.8
$\pm $ 10.5 mm for v-cut), which is also observable in the trajectories (see Fig. 3B). Foot markers, such as the ToeL, which is placed at the head of the fifth metatarsal, were tracked more closely by the inverse methods than by the tracking simulations. This can for example be seen for the right foot in the screenshot of sample 31 in Fig. 3A or in the trajectories of ToeL in Fig. 3B. For the GRFs, marker tracking showed slightly smaller RMSDs compared to coordinate tracking (*e.g*., 3.9
$\pm $ 2.6 BW% *vs*. 5.8
$\pm $ 2.6 BW% for v-cut). While there was no clear trend for the translation, coordinate tracking was generally closer to the reference angles, which were obtained with inverse kinematics, than marker tracking (*e.g*., RMSDs of 1.4
$\pm $ 1.8 deg *vs*. 3.5
$\pm $ 4.0 deg for v-cut). Especially the joint angles of the metatarsophalangeal (mtp) joint deviated more from the reference for marker tracking but also for coordinate tracking (see Fig. 3B). Inverse methods resulted in higher plantarflexion of the mtp joint compared to the tracking simulations, and coordinate tracking resulted in very high dorsiflexion during push-off. Furthermore, joint angles reconstructed with both tracking simulations were smoother than those obtained from inverse methods (see Fig. 3B). In principle, joint moments followed the same course for all reconstruction methods but showed different oscillations. The mean RMSDs of the joint moments were similar for the two tracking simulations but marginally smaller for marker tracking (*e.g*., 0.7
$\pm $ 0.7 BW BH% *vs*. 1.0
$\pm $ 1.1 BW BH% for v-cut). For both methods, arm joints showed smaller RMSDs for the moments than leg joints. The v-cut generally led to higher mean RMSDs compared to straight and curved running (see Table 3).

Motion | Inverse methods | Marker tracking | Coordinate tracking | Reference | |
---|---|---|---|---|---|

Marker in mm | SR | 6.7 $\pm $ 4.7 | 7.6 $\pm $ 4.8 | 13.4 $\pm $ 8.8 | Measured data |

CR | 6.8 $\pm $ 4.7 | 8.2 $\pm $ 5.5 | 13.5 $\pm $ 8.6 | ||

VC | 7.6 $\pm $ 5.1 | 9.5 $\pm $ 6.2 | 16.8 $\pm $ 10.5 | ||

GRF in BW% | SR | 3.4 $\pm $ 1.8 | 4.4 $\pm $ 2.4 | ||

CR | 4.2 $\pm $ 3.4 | 4.5 $\pm $ 2.3 | |||

VC | 3.9 $\pm $ 2.6 | 5.8 $\pm $ 2.6 | |||

Translation in mm | SR | 2.6 $\pm $ 1.1 | 2.7 $\pm $ 1.3 | Inverse methods | |

CR | 3.0 $\pm $ 1.7 | 2.9 $\pm $ 1.6 | |||

VC | 4.5 $\pm $ 1.9 | 3.2 $\pm $ 1.2 | |||

Angle in deg | SR | 2.3 $\pm $ 2.9 | 0.9 $\pm $ 0.7 | ||

CR | 2.5 $\pm $ 3.0 | 0.9 $\pm $ 0.7 | |||

VC | 3.5 $\pm $ 4.0 | 1.4 $\pm $ 1.8 | |||

Moment in BW BH% | SR | 0.6 $\pm $ 0.6 | 0.8 $\pm $ 0.8 | ||

CR | 0.7 $\pm $ 0.6 | 0.8 $\pm $ 0.9 | |||

VC | 0.7 $\pm $ 0.7 | 1.0 $\pm $ 1.1 |

Residual forces and moment of inverse methods and tracking simulations differed considerably. For inverse methods, the mean RMS residual pelvis forces and moments over all trials were 5.9 $\pm $ 2.3 GRF ${}_{\text{max}}$% and 0.9 $\pm $ 0.4 GRF ${}_{\text{max}}$ BH%, respectively. In contrast, marker and coordinate tracking simulations had maximum RMS residual pelvis forces of $6.0\cdot {10}^{-8}$ GRF ${}_{\text{max}}$%, pelvis moments of $1.7\cdot {10}^{-8}$ GRF ${}_{\text{max}}$ BH%, and joint moments of $1.5\cdot {10}^{-6}$ ${\mathrm{G}\mathrm{R}\mathrm{F}}_{\text{max}}$ BH%.

The scaled models, the experimental data, the result of the inverse methods, and the simulation results are provided online (Nitschke et al., 2022).

## Discussion

In this article, we showed that it is feasible to reconstruct measured motions by directly tracking marker and GRF data without task constraint in an optimal control simulation of a 3D full-body musculoskeletal model. We successfully tracked COD running motions without prior knowledge of the task nor the initial state. The presented formulation of the dynamics is therefore suited to reconstruct arbitrary motions. Marker tracking was superior to coordinate tracking and comparable to inverse methods in terms of marker errors while resulting in mutually and dynamically consistent kinematics and kinetics.

Marker tracking simulations took approximately three to four times longer to solve than coordinate tracking simulations (see Table 2) and therefore were computationally much more expensive. At the same time, marker tracking did not require considerably more CPU time per iteration, implying that the higher number of iterations mainly caused the large increase in computation time. The marker tracking has a higher complexity due to a higher non-linearity in the objective and gradients since marker positions must be obtained from the model states. In contrast, the generalized coordinates are part of the states and, therefore, optimization variables. Consequently, alternative optimization algorithms or formulations of the problem would have to be investigated to decrease the number of iterations and thus the computation time of marker tracking. In conjunction with this, the ratio of CPU time spent in the function evaluation of the NLP was higher for marker tracking since marker positions must be computed in every iteration.

Despite the recent advances, solving optimal control problems is computationally demanding, while inverse methods can be computed in seconds or minutes. However, the time required to solve an optimal control simulation highly depends on the formulation of the problem and its implementation. In this work, we used a large number of collocation nodes of up to 180 by reconstructing the motion at the original sampling frequency of 175 Hz, resulting in a large number of optimization variables and thus unknowns of up to about 80,000. In comparison, Venne et al. (2022) generated simulations with up to 106 nodes and 12,444 variables using multiple shooting. Furthermore, the choice of initial guess affects the time required for the optimization. We initiated the optimal control problem using a standing simulation to have an unbiased initial guess. Instead, the solution from inverse methods could be used, which decreases the time required since this initial guess is closer to the final solution. We chose not to do this because this would bias the simulations towards the result of the inverse methods and not allow for an independent comparison between the three methods. Therefore, optimization could be further accelerated by reducing the number of collocation nodes and using an informed initial guess depending on the application’s specific requirements.

Measured marker positions were tracked closest by inverse kinematics (see Table 3) since it estimates the kinematics for each time point separately without accounting for model dynamics. Neglecting dynamics makes inverse kinematics prone to track measurement noise and soft tissue artifacts, leading to smaller marker errors but also to high-frequency components in the movement (see Fig. 3B). Those high-frequency components make it necessary to filter estimated coordinates before using them as input for inverse dynamics. However, the choice of the cut-off frequency can considerably influence the resulting joint moments (Derrick et al., 2020). In contrast to inverse methods, optimal control simulation acts as a physical filter by accounting for model dynamics and minimizing effort in the objective. This eliminates the need to filter the data in advance but requires to balance tracking and effort in the objective to find a trade-off between close tracking and realistic neural excitation patterns.

Inverse kinematics tracked markers at the feet more closely than the simulations, even though marker positions were strongly affected by shoe deformation. The close tracking of the deformed and thus inaccurate marker positions resulted in unrealistically high plantarflexion of the mtp joint for inverse kinematics (see Fig. 3B). Similar to soft tissue artifacts, the deformations of the shoes are not modeled by the musculoskeletal model since virtual markers are rigidly attached to the model. However, in the optimal control simulations, high plantarflexion was prevented by modeling passive moments for all joints which became active in the mtp joint when plantarflexion exceeded 8 degrees (see supplemental information of Nitschke et al. (2020)). Consequently, the feet marker are then not tracked strictly, which results in higher marker errors for the simulation. In inverse kinematics and marker tracking simulation, the mtp joint angle could become more realistic by weighting the markers on the deformed shoe less than the other markers. In any case, weights should be adjusted only if it is appropriate considering the data, application, and biomechanical variables of interest since it might worsen the estimation of other variables and requires hand-tuning.

For marker and coordinate tracking simulation, input variables were tracked more closely than the variables not used in the objective of the optimization (see Table 3). Coordinate tracking follows the recorded marker data worst since errors propagate by tracking inaccuracies in the coordinates estimated with inverse kinematics. In detail, inaccuracies resulting from measurement, inverse kinematics, and coordinate tracking add up in contrast to marker tracking, where only the inaccuracies resulting from the measurement and tracking add up. Furthermore, errors made in the tracking of individual joint angles accumulate along the kinematic chain and result in a larger difference in the position of distal segments and thus of the marker positions. Therefore, distal segments like the hands and feet deviate considerably for coordinate tracking from the other two reconstruction methods (see Fig. 3A). These findings are in agreement with previous work for skeletal models (Febrer-Nafría et al., 2020; Venne et al., 2022) and proof that marker tracking is more accurate than coordinate tracking in terms of marker error. However, inverse kinematics which we used for comparison is not a gold standard since it is not reflecting the bone motion, but is subject to errors. Hence, estimated kinematics should be evaluated with bone pins or medical imaging. Regardless of the reconstruction method chosen, we strongly recommend analyzing the marker error carefully since this is the only available error measure when performing optical motion capturing.

Marker tracking reconstructed the measured GRFs slightly better than coordinate tracking (see Table 3). Even though we used the same weight of the GRF tracking term in both types of simulations (see Table 1), it might be that the GRF term had a higher influence on the overall objective in marker tracking than in coordinate tracking as the kinematic tracking data differed. Therefore, adjusting the weighting in the objective could counteract the different accuracies with respect to the GRFs, but would change the relation between GRF tracking and effort term.

The difference in GRF tracking between the marker and coordinate tracking might have also caused the slight difference in the reconstruction of joint moments. Again, it is necessary to note that the reference joint moments obtained with inverse dynamics do not represent a gold standard. However, joint moments can only be measured using instrumented implants. Alternatively to the comparison with inverse methods, the simulation could be evaluated by reconstructing simulated motions. Simulated data has the advantage that the ground truth would be known. Nevertheless, an evaluation with simulated data could hardly reflect all characteristics of real-world data like noise, errors, and soft tissue artifacts perfectly.

Although the constraint we introduced to remove the need for a task constraint (see Eq. (8)) allows the reconstruction of various motions, one minor downside is that it requires additional samples before the motion of interest. However, we also recommend analyzing longer time periods when applying inverse methods to reduce filtering artifacts. When reconstructing longer motions with simulation, it could be investigated to use moving horizon estimation (Bailly et al., 2021) where a new simulation is initiated using the end of the last simulation.

The analysis of residual forces and moments for inverse methods and tracking simulation highlights the difference between the two methods with respect to dynamic inconsistencies. In this study, inverse methods led to slightly higher residuals than recommended by Hicks et al. (2015). They recommend RMS residual forces lower than 5 GRF
${}_{\text{max}}$% and residual moments lower than 1% of GRF
${}_{\text{max}}$ times center of mass (CoM). Residuals in inverse methods could be reduced by manually adjusting the inertial parameters of the scaled model (Hicks et al., 2015). However, a manual adjustment to every participant is hardly feasible for large studies or in automated analysis pipelines. In contrast to inverse methods, the optimal control simulations had negligibly small residuals. As a result, estimated biomechanical variables are dynamically consistent. Therefore, there are no inconsistencies between the cause, *i.e*., neural excitation of the muscles, and the effect, *i.e*., the resulting motion. However, it has to be kept in mind that the simulation, in return, assumes that the dynamic model is perfect.

Potential inaccuracies and simplifications in musculoskeletal models can limit all reconstruction methods, *i.e*., inverse methods and optimal control simulation. In this study, we scaled the musculoskeletal model using marker positions of a static trial but did not personalize any muscle properties. For both simulation methods, it might be possible to further reduce tracking errors and improve muscle variable estimation when a better estimate of muscle parameters is available, for example, from strength tests (Hegarty et al., 2019) or medical imaging (Valente et al., 2017).

In the future, tracking marker positions directly instead of estimated generalized coordinates could offer further possibilities. It would, for example, allow personalizing model parameters based on measurement data within the optimal control problem by adding certain parameters (*e.g*., segment lengths) to the optimization variables instead of predefining them. Furthermore, simulations can be created even when marker data is incomplete, for example, due to occlusions (Venne et al., 2022). The time periods where marker data is missing can be excluded from the tracking objective since model dynamics and effort minimization will still produce a realistic movement for these periods. But most importantly, marker tracking simulations could also be driven by virtual marker positions extracted from video data, depth images, or radar technology instead of using marker-based optical motion capturing.

## Conclusions

In conclusion, we proved that it is feasible to directly drive optimal control simulations by marker and GRF data for 3D full-body musculoskeletal models to reconstruct COD running motions without estimating generalized coordinates in an intermediate step. We presented a detailed comparison of marker tracking simulations, coordinate tracking simulations, and inverse methods. In contrast to inverse kinematics and dynamics, optimal control simulation returns kinematics and kinetics, which are mutually and dynamically consistent. Dynamic consistency is especially important for the analysis of fast motions for example in sport science. Our results confirmed that marker tracking reconstructs measured marker positions more accurately than coordinate tracking. We, therefore, recommend using marker tracking simulations over coordinate tracking for reconstructive simulations, especially for applications investigating small changes in kinematics or kinetics. Nevertheless, coordinate tracking might still be advantageous when reference data is included in predictive simulations.

## Supplemental Information

### Frontal view of kinematics of a v-cut.

The result of inverse kinematics, coordinate tracking, and marker tracking is represented in red, orange, and green, respectively. The measured marker positions are displayed in blue. The motion had in total 149 samples at 175 Hz.

### Lateral view of kinematics of a v-cut.

The result of inverse kinematics, coordinate tracking, and marker tracking is represented in red, orange, and green, respectively. The measured marker positions are displayed in blue. The motion had in total 149 samples at 175 Hz.