Combined Effects of External Moments and Muscle Activations on ACL Loading during Numerical Simulations of a Female Model in OpenSim

Abstract

Background: Anterior cruciate ligament (ACL) injuries have been studied using a variety of methods and tools. However, each is hindered by specific limitations with respect to its application. Aim: To assess the combined effects of external moments and muscle activations on ACL loading using serial, forward dynamics (FD) simulations of single leg, hyperextension landings in OpenSim. Methods: The FD tool of OpenSim was iteratively run using different combinations of knee-spanning muscle activation levels, internal rotation and valgus knee moment magnitudes. A regression was conducted on the data in order to predict ACL loading under different conditions. Results: A purely abduction moment leads to greater mean ACL loading than a purely internal rotation moment or any combination of the two. Additionally, the generalized boosted regression model using both external moments and certain knee muscles identified the internal rotation moment as the most important variable in predicting the ACL load (R² = 0.9; p < 0.0001). Conclusion: This study demonstrated a novel and practical application of an OpenSim musculoskeletal model that supports the ACL injury mechanism of landing with low knee flexion angles, high muscle forces of the Quadriceps muscles and an external knee valgus moment, though further investigation is needed.

1. Introduction

Ruptures of the anterior cruciate ligament (ACL), one of the most serious injuries in sports, affect female athletes at an incidence rate 2.2 times higher than males [1]. Prevention of ACL injuries is the ideal goal as fewer than half of competitive athletes return to the same or higher activity level two years post-surgery following an ACL rupture [2]. ACL tears often occur in a non-contact manner when an individual quickly decelerates to change direction or lands from a jump [3]. About two thirds of ACL injuries occur during a knee valgus or varus injury mechanism [3], and the mechanics of the landing and cutting maneuvers that may result in these injuries are thoroughly studied [4,5,6]. However, almost one-third of the injuries occur during a hyperextension event [3], and this kind of injury has hardly been studied at all.

The forces that tear the ACL are well understood from cadaveric studies that apply various loads to knee specimens in order to assess their effects on ACL strain [7]. Such studies are valuable as measurements taken from strain transducers implanted in the ACL allow for the effects of potential loading mechanisms to be studied directly. As a result, the multi-planar loading of the knee with abduction and internal rotation moments has been shown to result in higher ACL strains compared to that of uni-planar loading [8]. These findings have been used in motion capture studies to identify potential risk factors for ACL injury [5,6,9]. However, both cadaver and motion capture studies have inherent limitations. Cadaver studies are limited in the number of tests that can be conducted due to wear on the specimens and the struggle to account for the role of the muscles during dynamic loading. Motion capture studies are limited to studying movements where no ACL injury occurs, and the forces acting on the ACL are largely unknown. Computer simulation studies (in silico) are well suited to exploring the complex interactions of external moments and muscle forces that produce high ACL strain values during injurious situations such as a single leg hyperextension landing.

Through programs such as OpenSim, forward dynamic (FD) simulations can be used to study various aspects of gait under conditions not typically available in experimental approaches [10]. Model parameters such as ligament properties and muscle activations can be manipulated at will, and the number of tests that can be performed is limited only by computation time. As a result, the interest in musculoskeletal modeling has increased in recent years, with researchers advancing towards anatomically correct knee models [11,12,13].

Although previous studies have used OpenSim modeling to assess knee kinematics under various conditions [14,15,16], no study has used these models for dynamic serial simulations to assess multi-planar loading during knee hyperextension on the ACL despite having been identified as a potential injury mechanism and a predictive tool for ACL injuries [17]. Therefore, the aim of this study was to explore the effects of various external moments and muscle activations on ACL loading using serial FD simulations of single leg, hyperextension landings in OpenSim.

2. Materials and Methods

2.1. Model

The model used in this research was a derivation of the Schmitz model (runner18_wholebody), which is available for use on the OpenSim web page [18]. This model represents a 77.5 kg female and includes a discrete element knee with a 6 degree of freedom tibiofemoral joint and 18 ligament bundles in the right leg, including 2 that represent the ACL [12].

The default initial pose of the model throughout all simulations was an upright, standing position with 0 degrees of hip and ankle flexion and 25 degrees of knee flexion, which has been used in previous cadaveric model studies [19,20,21]. All model coordinates were locked during the simulations except for the flexion/extension, internal/external rotation, abduction/adduction, and translations (Tx, Ty, Tz) of the right knee. In order to prevent nonphysiological movement in these directions, Coordinate Limit Force actuators were added to the model to constrain their ranges of motion (

Table 1. Default value and range for 6 DOF knee model.

Knee Coordinate	Default Value	Range
Angle	−25°	−110 to 0°
Adduction	0°	−7 to 7°
Rotation	0°	−25 to 25°
Tx	−0.006 m	−0.024 to 0.012 m
Ty	−0.449 m	−0.449 to −0.444 m
Tz	0 m	−0.003 to 0.003 m

). This type of actuator is inactive until a coordinate exceeds a user-defined range, at which point it supplies an inhibitory force to limit further movement outside of that range. Because the coordinate governing the patellofemoral joint was locked during the simulations, the attachment points of the quadriceps muscles were extended to include the path of the patella tendon in order to maintain their biomechanical function as knee extensors.

2.2. Force Conditions

External forces and moments were applied to the model to simulate dynamic loading conditions. A constant force of 1300 N in the upward (+Y) direction acting on the center of mass of the right tibia was added to simulate contact with the ground and was approximately two times body weight [22]. Knee internal rotation and valgus moments, whose maximum values were 90 N*m and 150 N*m, respectively, were varied in magnitude and applied to the model in different combinations, similar to previous cadaver and simulation studies [8,22]. Muscle activations of the hamstrings and quadriceps were also varied and tested in different combinations. In order to reduce the computation time and computer memory needed to perform the simulations, constant values for the moments and muscle activations were used throughout the entirety of each simulation. A list of forces applied to the model and their values are presented in

Table 2. External forces (a) present during musculoskeletal simulations and muscle parameters (b) governing the model.

a. Model Parameter	Value Range
Ground reaction force	1300 N
Rotation Moment	[0, 29.7, 59.4, 90] N*m
Adduction Moment	[0, −9.5, −99, −150] N*m
Muscle Activations	[0.2, 0.5, 0.8]
b. Muscle	Maximum Isometric Force (N)	Activation Range
Semimembranosus	1288	[0.2, 0.5, 0.8]
Semitendinosus	410	[0.2, 0.5, 0.8]
Biceps Femoris–long head	896	[0.2, 0.5, 0.8]
Biceps Femoris–short head	804	[0.2, 0.5, 0.8]
Rectus Femoris	1169	[0.2, 0.5, 0.8]
Vastus Medialis	1294
Vastus Intermedius	1365
Vastus Lateralis	1871
Sartorius	156	0.2
Gracilis	162
Medial Gastrocnemius	1558
Lateral Gastrocnemius	683

.

2.3. Forward Simulation

A custom MATLAB (MathWorks Inc., Natick, MA, USA) script based on methods by Mansouri and Reinbolt (2012) [23] was written to iteratively run the FD tool of OpenSim 4.0 Table 2a The script begins by defining the sample space for the parameter combinations of the activation levels for the quadriceps, hamstrings, and knee actuators. For each combination, the controls file is updated with values for each respective variable, and the FD tool is initialized. In OpenSim, the forward dynamics tool solves a set of differential equations that model musculoskeletal dynamics. The results of the integration describe how the states of the model (i.e., coordinates, coordinate velocities, muscle activation, muscle fiber length) change over time based on the forces and moments acting upon it (i.e., controls) [10]. Equilibrium actuator states are solved, and a force reporter analysis is used to obtain either the scalar force or torque for actuators or the calculated tension along the fiber path for muscles and ligaments for the duration of the simulation. From these results, kinematic/kinetic data of interest are gathered and recorded, and the next parameter set is tested. Based on the previous literature and the scope of this study, the following variables were recorded: maximum ACL loading from each bundle (anterior and posterior) and knee coordinates (Table 1) and knee muscle outputs (Table 2b) at the moment of maximum ACL loading for both bundles.

Input files for the simulations included the model file, a force set file, which contained the parameters of all point and torque actuators applied to the model, and a controls file, which defined both the duration of the simulated movement and the activation levels of each muscle and actuator during that simulation. Video analysis has shown that ACL injuries are likely to occur during the first 40 ms after initial contact [24], so a simulation time of 100 ms was deemed sufficient. In total, 3888 simulations were run. The data was stored in a Microsoft Office Excel spreadsheet for further analysis.

2.4. Statistical Analysis

From each simulation condition, the peak ACL force for each bundle was extracted. The peak muscle force before the peak ACL force was extracted. A temporal relationship was assumed where it would take an unknown amount of time from muscle contraction until the force was fully transferred to the ACL. Kinematics at the time of peak ACL force were extracted.

The inputs to the OpenSim model were the forces applied to the knee via moments and muscle forces. The forces in turn influence kinematics, which dictate the ACL load. Generalized linear models [25] were constructed to reflect this causality. The first model was the influence of kinematics on ACL forces. The second model was the influence of external moments on ACL forces. The third model combined external moments and muscle forces.

Because the relationship between forces and ACL loads may be non-linear and rely on complex interactions, a generalized boosted regression (GBM) was also created [26]. GBM is a machine learning algorithm that uses Friedman’s Gradient Boost algorithm to iteratively build a regression model via decision trees [27]. The algorithm begins by initializing a constant value upon which the model will be built. Then, residuals are calculated using the constant as the predicted value, and regression trees are made to fit those residuals. The model is updated by adding the constant value to the new regression tree multiplied by a learning rate. This process of making new regression trees based on the residuals of previous models and updating that model is repeated for a user-defined number of iterations or until additional regression trees fail to improve the model. The number of decision trees used to train the model in this study was set to 350, and the size of each tree could have a depth no larger than three. To prevent overfitting, cross-validation was implemented in the model with only a fraction (70%) of the data used to create each new tree in the model.

Fit between models was compared using the squared correlation coefficient (R2), and variable importance was compared using the adjusted partial correlation coefficient. As calculations of partial correlations are not possible for a GBM model, variable importance is calculated as the loss of root mean square when the variable is excluded. Due to the exploratory nature and large amounts of data available, alpha was set at 0.00001, and no p-value adjustments were performed.

3. Results

Descriptive statistics for the model outputs are shown in

Table 3. Mean and standard deviation of the ACL load for the anterior (a) and posterior (p) ACL bundles during each combination of external knee moments. Corresponding means and standard deviations of knee extension (Ext), abduction (Abd), rotation (Rot) and anterior translation (TX) are also recorded for each bundle.

		Mean & Std for aACL										Mean & Std for pACL
Rot Moment (N*m)	Add Moment (N*m)	aLoad (N)		aAngle (°)		aAbd (°)		aRot (°)		aTX (m)		pLoad (N)		pAngle (°)		pAbd (°)		pRot (°)		pTX (m)		Total Load (N)
0	0	527.1	72.1	−8.9	14.8	0	0.7	4.5	4.6	0.013	0.003	1129.4	190.6	−3.8	11.6	0.1	0.9	4.7	4.4	0.012	0.003	1656.5	262.7
29.7	0	590.2	74.1	−0.7	10.4	−0.6	0.3	9.5	3.2	0.012	0.002	1272.1	156.4	2.2	3.4	−0.7	0.2	10.3	1.7	0.011	0	1862.4	230.4
59.4	0	486.8	131.9	−7	15.4	−1.7	1.2	15.5	8.7	0.012	0.004	1087.8	231.1	3.6	2.2	−2.8	0.8	22.7	3.9	0.008	0.003	1574.6	363
90	0	369.9	65.9	−5.2	13.8	−2.7	1.3	21.2	9	0.01	0.004	921.1	118.1	3.6	1.5	−3.5	0.3	26.5	1.1	0.007	0.002	1291	184
0	−49.5	571.7	72.1	0.1	7.5	−0.8	0.3	7.8	2.5	0.012	0.001	1245.8	149.5	0.8	5.5	−0.8	0.2	7.9	2.2	0.011	0.001	1817.6	221.7
29.7	−49.5	585.7	105.9	3.2	1.8	−1.7	0.9	12.6	2.3	0.01	0.001	1274.7	175.8	3.2	1.8	−1.7	0.9	12.6	2.3	0.01	0.001	1860.5	281.7
59.4	−49.5	483.4	117.4	3.9	0.8	−3.7	0.2	24.2	1.5	0.008	0.002	1095.6	182.9	4	0.8	−3.7	0.2	24.2	1.5	0.007	0.002	1579	300.3
90	−49.5	432.2	95.9	3.7	1.1	−3.9	0.3	26.2	1.4	0.007	0.002	1017.2	149.5	3.9	0.9	−3.9	0.3	26.2	1.4	0.007	0.002	1449.4	245.5
0	−99	597.5	60.4	1.9	4.2	−1.1	0.2	8.6	1.6	0.011	0.001	1294.8	110.7	2.2	2.6	−1.1	0.2	8.7	1.4	0.011	0	1892.3	171.2
29.7	−99	618.2	108	3.3	4.3	−3.1	0.6	13	3.5	0.01	0.002	1316.4	177.5	3.8	1.3	−3.2	0.6	13.3	3.5	0.01	0.002	1934.6	285.5
59.4	−99	481.4	76.6	4.2	0.9	−4.2	0.2	23.4	2.1	0.007	0.001	1093.6	123.2	4.3	0.8	−4.2	0.2	23.5	2.1	0.007	0.001	1575	199.8
90	−99	473.4	77.2	4.2	0.7	−4.5	0.2	26.1	1	0.007	0.002	1076.1	120.1	4.2	0.7	−4.5	0.2	26.1	1	0.007	0.002	1549.5	197.3
0	−150	667	99.1	1.6	5.4	−3.6	1	8.4	2.2	0.012	0.001	1382.1	168.5	2.3	1.9	−3.7	0.8	8.8	1.8	0.012	0.001	2049.1	267.6
29.7	−150	666.2	96.9	3.2	3.7	−5.3	1.6	12.4	3.5	0.011	0.002	1374	155.7	3.5	1.1	−5.3	1.5	12.4	3.2	0.011	0.001	2040.2	252.6
59.4	−150	558.5	64.5	3.3	6.3	−6.1	1.3	20.7	3.3	0.008	0.002	1200.2	101.3	4.5	0.7	−6.2	1	21.4	1.7	0.008	0.001	1758.7	165.8
90	−150	499.7	38	−0.6	9.1	−6.1	1.8	24.9	4.6	0.01	0.003	1108.4	69.5	4.6	1.6	−6.9	1.2	25.6	1.1	0.008	0.002	1608.2	107.5

. Positive values for kinematics represent knee extension, adduction, internal rotation, and anterior translation. The combination of 0 N*m of internal rotation moment and 150 N*m of valgus moment resulted in the highest mean loading for both bundles (a: 667.0 ± 99.1 N, p: 1382.1 ± 168.5 N), while the combination of 90 N*m of rotation and 0 N*m of valgus resulted in the lowest mean loading for both bundles (a: 369.9 ± 65.9 N, p: 921.1 ± 118.1 N).

The ACL was modeled as two bundles in OpenSim: the anterior and posterior bundles. There was no difference in magnitude or direction between each independent variable and the load on each respective bundle. For simplicity, results of the regression analysis are therefore presented only for the anterior bundle.

The results of the generalized linear regression using only knee kinematics are presented in

Table 4. Regression table for a generalized linear model using knee kinematics to predict ACL loading.

Variable	Coefficient	Partial Correlation	p-Value
Flexion	28.92	0.78	<2 × 10−16
Rotation	−4.38	−0.47	<2 × 10−16
Valgus	4.3	0.11	<2 × 10−11
Anterior Translation	46,300	0.86	<2 × 10−16
Medial Translation	−58,880	−0.53	<2 × 10−16

. The model R2 was 0.91 and of the kinematic variables the anterior tibial translation had the greatest partial correlation (0.86).

The results of the generalized linear regression using only external moments are presented in

Table 5. Regression table for a generalized linear model using externally applied knee moments to predict ACL loading.

Variable	Coefficient	Partial Correlation	p-Value
Int. Rotation Moment	−1.85	−0.55	<2 × 10−16
Valgus Moment	0.67	0.37	<2 × 10−16

. The model R2 was 0.37, and the rotation moment had a greater partial correlation than the valgus moment.

The results of the generalized linear regression for all forces are presented in

Table 6. Regression table for a generalized linear model using hamstring and quadriceps muscles and externally applied knee moments to predict ACL loading. The final column of the table was calculated using the same predictor variables in a generalized boosted model (GBM).

Variable	Coefficient	Partial Correlation	p-Value	Root Mean Square Loss
Valgus Moment	0.82	0.41	<2 × 10−16	55
Int. Rotation Moment	−3.16	−0.64	<2 × 10−16	122
Biceps Femoris Long	−0.05	−0.15	<2 × 10−16	51
Biceps Femoris Short	−0.01	−0.03	<2 × 10−16	42
Rectus femoris	1.38	0.13	<2 × 10−16	57
Semimembranosus	−0.07	−0.25	<2 × 10−16	53
Semitendinosus	−0.11	−0.18	<2 × 10−16	46
Vastus intermedius	0.68	0.13	<2 × 10−16	47
Vastus Lateralis	−1.26	−0.3	<2 × 10−16	56
Vastus Medialis	0.25	0.04	<2 × 10−16	63

. The model R2 was 0.59, and the muscle with the highest partial correlation was the Vastus Lateralis muscle at −0.30.

The results of the GBM are presented in Figure 1 and the last column of Table 6. The model R2 was 0.9. The most important variable (Table 6, last column) was the internal rotation moment followed by the Vastus Medialis muscle.

4. Discussion

The aim of this study was to explore the effects of external moments and muscle activations on ACL loading using serial FD simulations of single leg, hyperextension landings in OpenSim. The main finding of the study was that external knee moments alone explained much of the variability in ACL forces during single leg landings. Furthermore, the inclusion of muscles in the model increased the model fit by 0.22.

The valgus moment was positively correlated with ACL loading, while the internal rotation moment was negatively correlated with ACL loading. This finding somewhat contradicts the results of Kiapour et al., who showed that multiplanar loading of cadaver knees produced the highest strain on the ACL [8]. These differences may be due to the effects of joint geometry, which was modeled in our study as a flat surface controlled by coordinate limit force actuators. Similar to our findings, Kiapour et al. found that in a regression, the valgus moment explained more variability in ACL strain than the internal rotation moment [8]. A cadaveric study has found that an internal rotation moment increases ACL forces under hyperextension [28] but without weight bearing or muscle forces included. A post-hoc exploration of the data revealed that the internal rotation moment correlated strongly with all four heads of the Quadriceps muscles. Complex interactions may explain why cadaveric models that produce clinically relevant ACL injury patterns have only recently been created [29]. The high R2 achieved with the BGM model supports the importance of complex interactions between external moments and muscle forces.

The risk of ACL injury increases with decreased knee flexion upon landing [5]. The results of our study support this finding, as the knee extension angle had a high partial correlation with ACL forces (Table 3), second only to anterior tibial translation. Combined with the high activation of the quadriceps muscles and the reduced ability of the hamstrings to counteract anterior tibial translation due to the poor moment arm with small knee flexion angles, landing close to full knee extension can result in very high loads on the ACL. These findings are consistent with recent studies using biplane fluoroscopy during in-vivo movements that have found relationships between patellar tendon strain and ACL strain [30] as well as between greater knee extension angles and ACL strain [31].

The present study is not without its limitations. Because the Schmitz model, which includes 16 ligament bundles in addition to the 2 that represent the ACL, was used in this study, it is possible that other ligament bundles were being loaded more under certain conditions, thereby alleviating the loading on the ACL. When preparing knee specimens for cadaver studies related to the ACL, other knee ligaments are typically removed to nullify their influences on the results. Therefore, the Schmitz model was chosen because it is one of the most anatomically complete knee models with 6 degrees of freedom in order to create as realistic of a simulation as possible. Because the use of different models would likely result in different results [16], using the most accurate model is important. Furthermore, because the Schmitz model represents a female, the results of this study may only be generalizable to that population.

The inclusion of Coordinate Limit Force actuators and the use of constant values for the moment and muscle activations represent the tradeoff between computation time and biofidelity commonly found in musculoskeletal modeling. These decisions were made in order to reduce the computation time and computer memory needed perform the simulations. However, the authors acknowledge that they may also underscore the full effect of the interactions within the musculoskeletal system when exposed to external moments.

Lastly, in order to limit the focus of this study to only knee kinematics/kinetics, all other model coordinates, including the coordinate governing the patellofemoral joint, were locked during the simulations. As a result, musculotendinous structures attached to this joint, namely the patella tendon and the quadriceps muscles, would have had no contribution to knee kinematics/kinetics during the simulations. In order to maintain their biomechanical function as the knee extensor, the attachment points of the quadriceps muscles were extended to include the path of the patella tendon. While not anatomically accurate, this change represents a close approximation, which the authors believe minimally impacts the results of this study.

5. Conclusions

This study demonstrated a novel and practical application of an OpenSim musculoskeletal model to better understand knee kinetics/kinematics during human movement by mimicking the methodology of ACL cadaver studies. Our results support an ACL injury mechanism of landing with low knee flexion angles, high muscle forces of the quadriceps muscles and an external knee valgus moment. Furthermore, our results propose a protective effect of the external tibial internal rotation moment on ACL loading under these conditions. This result was unexpected and may warrant further investigation into its role in ACL loading.