An approach to generate noncontact ACL-injury prone situations on a computer using kinematic data of non-injury situations and Monte Carlo simulation

Abstract

ACL-injuries are one of the most common knee injuries in noncontact sports. Kinematic data of injury prone situations provide important information to study the underlying ACL-injury mechanisms. However, these data are rare. In this work an approach is presented to generate injury prone situations for noncontact ACL-injuries on a computer. The injury prone situations are generated by a musculoskeletal simulation model using kinematic data of a non-injury situation and the method of Monte Carlo simulation. The approach is successfully applied to generate injury prone landings in downhill ski racing. The characteristics of the obtained injury prone landings are consistent with video recordings of injury cases.

Keywords:

• ACL-injuries
• musculoskeletal simulation
• Monte Carlo simulation
• optimal control

Introduction

Anterior cruciate ligament (ACL) injuries are one of the most common knee injuries in sports (Ireland 1999; Griffin et al. 2006). About 70% of ACL injuries occur without physical contact between athletes. Typical situations, which result in noncontact ACL injuries, are awkward landings, sidestep cutting or pivoting (Boden et al. 2000; Griffin et al. 2006). ACL injuries imply high surgical costs, a long and intensive rehabilitation time and have a devastating influence on a patient’s activity level and quality of life (Yu and Garrett 2007). Therefore, the prevention of noncontact ACL injuries is of utmost importance. An important step regarding the development of effective injury prevention programs is a detailed understanding of the underlying injury mechanism and related risk factors (Bahr and Krosshaug 2005).

If video sequences of ACL injury cases are available, video analysis studies provide important information on the athlete’s kinematics (e.g., joint angles, orientation of the trunk) during the events leading to injuries and the actual ACL injury situations. Video analysis has been applied to study injury cases in different sports such as basketball, handball, football and skiing (Carlson et al. 2016). However, video sequences of actual injury situations are rare and video analysis studies are restricted to kinematic analyses.

Due to ethical reasons experimental studies, which replicate possible injury situations, are not feasible. However, human musculoskeletal simulation models offer the possibility to simulate injury prone maneuvers on a computer (e.g., Gerritsen et al. 1996; McLean et al. 2003; Lin et al. 2012; Heinrich et al. 2014; Weinhandl and O’Connor 2017). The simulation results provide information on variables that cannot be measured in an experiment or obtained by video analysis, such as joint moments, joint power, muscle activation as well as muscle and ligament forces (Pandy and Andriacchi 2010). In addition, potential risk factors can be investigated. Selected input variables to the simulation model can be varied and the corresponding effect on the loading of the knee joint can be determined.

Musculoskeletal simulation models typically follow a two-step approach to generate and subsequently analyse possible injury situations (e.g., Gerritsen et al. 1996; McLean et al. 2004, 2008; Heinrich et al. 2014; Weinhandl and O’Connor 2017). In a first step, the simulation model is applied to reproduce measurement data of a well-coordinated and safe situation and to validate the model. In a second step input variables are perturbed to simulate possible injury-prone situations. While input variables were manually perturbed in early simulation models (e.g., Gerritsen et al. 1996), recent models are combined with a Monte Carlo (MC) method to investigate injury-prone situations (McLean et al. 2004, 2008; Lin et al. 2012; Heinrich et al. 2014). The MC method is a numerical method to investigate the response of a system under random input. Thereby, the input variables and parameters are randomly varied n-times and the systems response is computed and evaluated statistically.

In the simulation studies of McLean et al. (2004, 2008) a musculoskeletal simulation model was combined with the MC method to generate repeated forward dynamic simulations with randomly perturbed initial states and perturbed muscle excitation patterns. This methodological approach is well suited to explore the global solution space and, consequently, worst case scenarios. However, this approach does not guarantee that the simulations result in physiologically reasonable locomotion because task-oriented adaptive behavior of the athlete is not considered. During an out-of-balance landing maneuver, for example, an athlete might aim for a recovery movement to regain balance. This behavior can be considered adding task constraints to the simulation and formulating a corresponding optimal control (OC) problem that minimizes a given objective function. The OC problem can be adapted to consider the task-oriented adaptive behavior of an athlete (Heinrich et al. 2014).

In this work we present a computational approach, which allows one to generate noncontact ACL-injury prone situations based on a musculoskeletal simulation model, kinematic data of a non-injury landing and the MC method, and considers task-oriented adaptive behavior of the athlete. The developed approach is applied to generate injury prone landings corresponding to the so-called landing back-weighted situation, which was identified as one of three main ACL injury situations in World-Cup alpine downhill ski racing (Bere et al. 2011).

Methods

A two-step approach is presented to generate injury-prone jump landing maneuvers in competitive alpine downhill skiing. In a first step, a reference simulation was computed tracking measured kinematic data of a skier performing a safe jump landing maneuver. In a second step, a series of perturbed jump landing maneuvers was generated and injury-prone landing maneuvers were identified.

Musculoskeletal model

The planar musculoskeletal multibody simulation model of an alpine skier with two flexible skis described thoroughly in Heinrich et al. (2014) was applied (Figure 1). The skier consisted of seven rigid segments and the flexible skis were treated as a multibody system each with nine rigid segments, which were connected with rotational springs and dampers. The resistance of the ski boot with respect to the ankle dorsiflexion and plantarflexion was modeled with a passive joint moment acting at the ankle joint (Eberle et al. 2017). The dynamics of the multibody ski-skier model was given by the following equation (1) $$f(q,\dot{q},\ddot{q},F_M(L_{CE},{\dot{L}}_{CE},a\left)\right)=0.$$(1)

Here q, $\dot{q}$ and $\ddot{q}$ denote the generalized coordinates, velocities and accelerations, and F_M the muscle forces, which depend on the muscle fiber lengths L_CE, the contraction velocities ${\dot{L}}_{CE}$ and muscle activations a. The movement of the skier was controlled by 16 muscles—eight muscles for the right and left lower extremity, respectively: iliopsoas, glutei, hamstrings, rectus, vasti, gastrocnemius, soleus, tibialis anterior (Figure 1). Each muscle was modeled with a three element Hill-type model (Zajac 1989) with contraction dynamics (McLean et al. 2003) (2) $$g(L_{CE},{\dot{L}}_{CE},a,q)=0$$(2) and activation dynamics (He et al. 1991) (3) $$h(a,\dot{a},u)=0$$(3) where u denotes the muscle excitations. The muscle properties were taken from Gerritsen et al. (1996). Combining equations (1)(1) $$f(q,\dot{q},\ddot{q},F_M(L_{CE},{\dot{L}}_{CE},a\left)\right)=0.$$(1) , (2)(2) $$g(L_{CE},{\dot{L}}_{CE},a,q)=0$$(2) and (3)(3) $$h(a,\dot{a},u)=0$$(3) yielded the dynamics of the musculoskeletal model as a system of implicit differential equations (van den Bogert et al. 2011) (4) $$I(x,\dot{x},u)=0$$(4) with 82 state variables $x=[q,\dot{q},L_{CE},a]$ (25 generalized coordinates q of the skier-ski model, 25 generalized velocities $\dot{q}$, 16 muscle fiber lengths L_CE and 16 muscle activations a) as well as 16 control variables u.

Figure 1 Musculoskeletal skier with two legs and flexible skis. Blue lines depict trunk and right leg, green lines left leg. Muscles are depicted for right extremity (in red colour).

Reference landing

In the present application example a reference landing was generated by a tracking simulation of a safe landing movement in downhill skiing (Nachbauer et al. 1996). For this an OC problem was formulated. The task of the OC problem (OPref) was to find state trajectories x(t) and control trajectories u(t) such that a given objective function J (5) $$J=J_1+J_2$$(5) was minimized, while subjected to the following constraints (6a) $$I\left(x\right(t),\dot{x}(t),u(t\left)\right)=0,\hspace{1em}0\le t\le T$$(6a)

(6b) $$0\le u\left(t\right)\le 1,\hspace{1em}0\le t\le T$$(6b) (6c) $$x_{hip}\left(t\right)=x_{hip,\mathit{\text{data}}}\left(t\right),\hspace{1em}t=0.$$(6c)

The constraints referred to the dynamics of the musculoskeletal skier-ski model (6a) and boundaries on the control variables u (6b), because muscle excitations are limited to unity bound constraints (Erdemir et al. 2007). In addition, the initial position of the hip of the skier $x_{hip}(0)$ was constrained to match the initial measured hip position $x_{hip,\mathit{\text{data}}}(0)$, such that the jump height of the skier was not altered in the simulation (6c).

The objective function J in the OC problem (OPref) consisted of two terms J₁ and J₂. Through the first term J₁ the N generalized coordinates of the musculoskeletal model $q_i(t),i=1,\dots ,N$ were fitted in a weighted least square sense to the data of the measured reference movement $q_{i,\mathit{\text{data}}}(t),i=1,\dots ,N$ scaled by factors $1/{\sigma }_i,i=1,\dots ,N$

(7) $$J_1=\frac{1}{N·T}{\int }_0^T\sum _{i=1}^N{\left(\frac{q_i\left(t\right)-q_{i,\mathit{\text{data}}}\left(t\right)}{{\sigma }_i}\right)}^{2}dt.$$(7)

The second term J₂ was used to resolve muscle redundancy (Laughlin et al. 2011; Mokhtarzadeh et al. 2013) where the term J₂ represented the sum of squared muscle activation patterns $a_i(t)$ averaged across the M muscles and the simulation time T and scaled by the weight factor w_mus

(8) $$J_2=\frac{w_{mus}}{M·T}{\int }_0^T\sum _{i=1}^M{a}_i{\left(t\right)}^{2}dt.$$(8)

The OC problem (OPref) was solved by the method of direct collocation (Betts 2010). The implicit midpoint rule was applied to discretize the system dynamics and to transform the OC problem (OPref) in a constrained Nonlinear Programming Problem (NLP) (van den Bogert et al. 2011). The NLP was solved efficiently by a primal-dual interior point optimization method (IPOPT, Wächter and Biegler 2006). For this analytical derivatives of the objective function and constraints were provided and the sparsity structure of the Jacobian of the constraints were exploited (van den Bogert et al. 2011).

Perturbed landing simulations

Based on the reference simulation the MC method was applied to generate a series of n = 1000 landing simulations, where the initial posture of the skier prior to ground contact was randomly perturbed. Regaining a balanced position represents the behavior that would be expected from a real skier who finds him/herself in a perturbed position prior to ground contact. Consequently, an OC problem was formulated to represent task-oriented adaptive behavior of the athlete and guarantee that the simulations result in physiologically reasonable locomotion.

Specifically, the posture of the skier was randomly perturbed 100 ms before touchdown (time t₁) by adding random numbers to the joint angles of both legs and trunk orientation. Touchdown was defined when the ground reaction force on the ski segment below the ski boot exceeded 10 N. Random numbers $({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_7)=ϵ$ were generated assuming Gaussian distributions with mean zero and standard deviations, which were derived from Barone et al. (1999) for the joint angles (hip: $9.4°$; knee: $8.7°$; ankle: $5.2°$) and the trunk orientation ($9.2°$).

With the initial posture of the skier $x_{pos}(t_1)$ constrained to the perturbed position in the reference simulation ($x_{pos}^{ref}(t_1)+ϵ$) an optimal control problem similar to the reference simulation was solved. In contrast to the reference simulation, the tracking term J₁ in the objective function, was replaced by a task constraint (10d) referring to the final position of the skier (at time T), where the skier was observed to be in a well-balanced downhill position. The reason for this was to generate predictive simulations where the skier attempts a recovery movement and aims for a well-balanced downhill position without constraining the skier’s movements during landing. In summary, the following OC problem (OPper) was formulated and solved in analogy to the reference simulation. Find state trajectories x(t) and control trajectories u(t) such that the objective function J (9) $$J=J_2$$(9) was minimized, while subjected to the following constraints (10a) $$I\left(x\right(t),\dot{x}(t),u(t\left)\right)=0,\hspace{1em}{t}_1\le t\le T$$(10a)

(10b) $$0\le u\left(t\right)\le 1,\hspace{1em}{t}_1\le t\le T$$(10b) (10c) $$x_{pos}\left(t\right)=x_{pos}^{ref}\left(t\right)+ϵ,\hspace{1em}t=t_1$$(10c) (10d) $$x_{pos}\left(t\right)=x_{pos}^{ref}\left(t\right),\hspace{1em}t=T.$$(10d)

Data analyses

For each landing simulation, ACL tensile forces were calculated using a data-driven knee model (Kernozek and Ragan 2008; Laughlin et al. 2011; Weinhandl and O’Connor 2017). The perturbed simulations with ACL tensile forces in the right or left knee higher than a specified threshold of 2000 N (McLean et al. 2008; Heinrich et al. 2014) were collected and classified as injury prone landings. The remaining perturbed simulations were classified as non-injury prone landings. The threshold of 2000 N was specified based on measured data of maximum tensile strength of ACL obtained by Woo et al. (1991).

The kinematics of injury prone landing maneuvers were compared quantitatively with non-injury prone landing maneuvers. A systematic video analysis of ACL injury cases during jump landings in downhill skiing showed that during the flight phase the skier lost balance backwards and landed asymmetrically on the ski tails with nearly extended knees (Bere et al. 2011). Consequently, the orientation of the trunk segment, the joint angles at the hip, knee, and ankle of the lower extremity corresponding to the peak force in the ACL as well as the differences in the joint angles at the hip, knee, and ankle of both lower extremities were compared. Instead of the joint angles of the second lower extremity, these differences were included in the analysis to account for asymmetries in the posture of the skier. Because the specified variables did not follow a Gaussian distribution, a Mann–Whitney U-test was used to investigate the differences of the kinematic posture variables between injury prone landings and non-injury prone landings. In addition, Spearman rank correlations were computed to examine the association between the kinematic posture variables and peak force in the ACL.

Results

In the reference simulation the skier landed in a slightly asymmetric position from a landing height of 1.29 m. Peak ACL force occurred 33 ms after initial ground contact and amounted to 919 and 293 N in the right and left knee, respectively. In the series of n = 1000 perturbed landing simulations 219 cases resulted in ACL tensile forces exceeding the specified threshold of 2000 N in the right knee and/or left knee (Figure 2 shows the distribution of peak ACL forces in the right and left leg corresponding to the 1000 perturbed landing simulations). Consequently, 219 cases were classified as injury prone landing maneuvers and 781 cases as non-injury prone landing maneuvers.

Figure 2 Distribution of peak ACL forces for right leg (blue bars) and left leg (yellow bars) in n = 1000 simulations.

At the time of initial ground contact, all specified kinematic posture variables except the difference in the ankle joint angle (ankle diff) showed a significant difference between the injury prone landings and the non-injury prone landing maneuvers (Table 1). Specifically, in the injury prone situations the trunk of the skier was more backward oriented, hip and knee flexion were decreased and ankle plantarflexion was increased. In addition, the landing position of the skier was observed to be more asymmetric, which was primarily caused by the difference in hip flexion between both legs (Table 1, Figure 3).

Figure 3 Box plots of the trunk orientation and joint angles of the leg corresponding to the peak ACL force (hip${}_{\mathrm{p}}$, knee${}_{\mathrm{p}}$, ankle${}_{\mathrm{p}}$) as well as the joint angles of the contralateral leg (subscript CL) in the injury prone cases (a) and the non-injury prone cases (b). The posture of the skier mimicking the corresponding median values is displayed in (c) and (d). For better visualization the snow surface is assumed to be horizontal and the posture of the skier is rotated backward by the slope inclination of 27.8°.

Table 1 Median as well as first and third quartile $[Q_1,Q_3]$ of the trunk and joint angles of the skier at initial ground contact (IGC) corresponding to the injury prone and non-injury prone situations (deg).

Variable	Injury prone	Non-injury prone	P-value
Trunk	−38.4 [−40.8, −35.5]	−50.6 [−57.3, −45.4]	<0.001
Hip	31.3 [27.9, 36.2]	43.8 [37.2, 50.1]	<0.001
Knee	9.7 [6.5, 13.7]	27.6 [19.4, 38.5]	<0.001
Ankle	8.9 [7.2, 11.6]	14.0 [12.9, 15.8]	<0.001
Hip diff	13.2 [7.0, 21.9]	0.2 [−7.9, 7.8]	<0.001
Knee diff	9.6 [3.8, 14.7]	7.1 [−2.0, 14.8]	0.003
Ankle diff	0.4 [−1.1, 1.7]	0.5 [−1.6, 2.5]	0.551

The variables Hip diff, Knee diff and Ankle diff were used to account for asymmetries in the posture of the skier.

The bivariate rank correlations showed significant relations between almost all kinematic posture variables and the peak ACL force (Table 2, p < 0.001). Only the differences in the ankle joint angles had a non-significant correlation coefficient (p = 0.148). In particular, increased backward lean, less hip, knee and ankle dorsiflexion as well as increased asymmetry in the hip and knee were associated with increased peak ACL force.

Table 2 Spearman rank correlations between the kinematics of the skier at initial ground contact (IGC) and peak ACL force.

Variable	r	P-value
Trunk	0.88	<0.001
Hip	−0.66	<0.001
Knee	−0.78	<0.001
Ankle	−0.68	<0.001
Hip diff	0.42	<0.001
Knee diff	0.20	<0.001
Ankle diff	0.05	0.148

The variables Hip diff, Knee diff and Ankle diff were used to account for asymmetries in the posture of the skier.

Discussion

In this work an approach was presented to generate ACL-injury prone situations by a musculoskeletal simulation model using kinematic data of a non-injury situation. In the approach the MC method was applied where the initial kinematic state of a reference situation was randomly perturbed and the system response was computed n times.

The presented approach was successfully demonstrated for jump landing maneuvers in downhill skiing. In the simulated injury prone landings, the skier got into a backward leaning position and landed asymmetrically on one ski tail with extended hips and knees (Table 1). The description of the obtained injury prone landings is consistent with other studies, where the same characteristics were observed (Bally et al. 1989; Bere et al. 2011; Heinrich et al. 2014).

The injury prone landings were generated based on an OC problem formulation combined with the MC method. In past studies focusing on a sidestep cutting maneuver (McLean et al. 2004, 2008) musculoskeletal simulation models combined with the MC method were applied to investigate certain injury situations. However, these studies involved forward dynamic simulations and neglected the task-oriented adaptive behavior of the athlete. The present formulation as an OC problem offered the possibility to consider the task-oriented adaptive behavior of the athlete and assure that the generated simulations result in physiologically reasonable locomotion. Especially, a task constraint on the final position of the skier was added. This formulation considered that an unbalanced athlete tried to regain balance and to recover the end position of the reference landing. Alternative neuromuscular control strategies or task constraints can be easily incorporated in the presented approach by changing the constraints or the objective function J or adding corresponding constraints to the OC problem, respectively.

The ordinary differential equations for musculoskeletal dynamics are often numerically stiff and highly nonlinear. Thus, the forward dynamic simulations require small time steps and OC problems are time-consuming to solve and have bad convergence properties. van den Bogert et al. (2011) presented an efficient method to solve forward dynamic musculoskeletal simulations and OC problems using an implicit formulation of the musculoskeletal dynamics. Although the OC problem (OP) can be solved by shooting or collocation methods, it has been shown that the direct collocation method is more efficient to solve OC problems in human movements (Porsa et al. 2016). Thus, a direct collocation method was applied to solve the OC problem (OPper), where we provided analytical derivatives of the objective function and constraints and exploited the sparsity structure of the Jacobian of the constraints (van den Bogert et al. 2011). Solving the OC problem (OPper) requires approximately 15 min on a common computer with an Intel® Core^TM i5-6600K CPU 3.50 GHz. Consequently, solving n OC problems in series requires approximately n times 15 min. The computation time of n simulations, however, can be reduced by computing the simulations in parallel on a single multi-core computer or a cluster of computers (parallelization) because the simulations of the MC method are independent of each other.

In this work m = 219 injury prone landings were collected in n = 1000 realizations. This implies that approximately every fifth jump is injury prone. The number of injury prone situations is affected by the reference landing, which is the basis of the presented approach. Schindelwig et al. (2015) investigated the landing height of jumps in ski racing. The landing height of the reference landing (1.29 m) was compared with similar jumps in Schindelwig et al. (2015) (MF, similar steepness and velocities) and is almost equal to the maximum landing height (1.37 m). In the OC problem formulation (OPper) the initial hip position was fixed such that the landing height was not altered. In Heinrich et al. (2014) the initial hip position was not fixed which resulted in a reference landing (baseline simulation) with a landing height of 0.9 m. Since the landing height of 0.9 m is closer to an average landing height (Schindelwig et al. 2015), the predicted number of injury prone situations was lower and more realistic. However, since we are interested in injury prone landings, a reference landing which is almost close to the maximal landing height is beneficial.

The presented approach has certain limitations. A limitation could be the objective function which was used to control the movements and resolve muscle redundancy. In particular, a control strategy was applied minimizing the sum of squared muscle activations. This control strategy is a commonly assumed muscle coordination strategy and has been used in recent simulation studies investigating jump landing movements (Laughlin et al. 2011; Mokhtarzadeh et al. 2013). To our knowledge control strategies in injury prone situations are not described in the literature. In the predicted injury prone situations the applied control strategy penalized co-contraction and produced high quadriceps muscle forces and low hamstrings co-contraction of the skier to regain balance. Interestingly, similar results were shown in Barone et al. (1999), who captured and analyzed an accidental ACL rupture during jump landing in skiing. There the muscle activity was measured by EMG measurements in an ACL-injury accident during a jump landing maneuver in alpine skiing. It was observed that the quadriceps muscle was activated prior to the impact whereas the hamstrings did not show relevant activity until the ski was completely clapped on the ski slope surface.

In conclusion, injury prone landing maneuvers in downhill skiing were successfully generated by the presented approach without any ethical justification problems. The computer generated injury prone landings can be used to study possible injury mechanisms, the effects of equipment (e.g. Eberle et al. 2017) or other risk factors (environmental, biomechanical, neuromuscular). The presented approach requires a musculoskeletal simulation model, data of a non-injury situation, an objective function and a specified threshold for the ACL-loads to generate kinematic and kinetic data by the MC method. Hence it might be suitable to generate simulations mimicking other typical noncontact ACL-injury prone situations such as sidestep cutting maneuvers, landing maneuvers in other sports, or stop-jump tasks by the presented approach. Based on the generated situations underlying injury mechanisms, risk factors or potential protective measures can be investigated on a computer. These results may provide important information for developing effective injury prevention programs aiming at the reduction of actual ACL injuries.

Disclosure statement

The authors state that there are no conflicts of interest.