Introduction

It is of great importance to know internal knee–ligament loading to develop better surgical procedures and rehabilitation regimens for anterior cruciate ligament (ACL)-deficient patients. Many researchers have used in vitro experiments to study the changes in anterior tibial translation (ATT) and ligament force distribution of the knee after transection of the ACL.2,14,22,25,26,42,43,46 They have shown increases in ATT and in situ forces in other soft-tissue structures, such as the medial collateral ligament (MCL) and the meniscus after removal of the ACL. Nevertheless, the external loads applied to the knee in these experiments were smaller than the loads during walking due to the limitation of the experimental systems, so these in vitro experiments could not replicate in vivo ligament loading during gait.

A few numerical models have been used instead to estimate in vivo ligament loading in the ACL-deficient knee during gait.27,40 Liu and Maitland27 developed a knee model in the sagittal plane to study knee mechanics of the ACL-deficient knee during level walking. This model was limited to a single selected position during early stance phase of gait. Shelburne et al. 40 developed a three-dimensional model of the lower extremity with a six degrees-of-freedom knee joint to calculate shear force and ligament loading in the ACL-deficient knee during walking, and compared the results to those of the healthy knee. In their study the joint angles, ground reaction forces, muscle forces, and joint-reaction forces of the ACL-deficient knee were assumed to be identical with those of the healthy knee. This assumption contradicted previous findings that ACL-deficient subjects have different gait patterns than healthy subjects.34 ACL-deficient knees have been reported to have different muscle activation strategies, such as strong co-contractions, and diminished vastus lateralis (VL) and lateral gastrocnemius (LG) muscle control.35,45 The previous numerical models did not estimate muscle forces based on measured muscle activations, i.e., electromyographic (EMG) signals, so they may have difficulty predicting the abnormal muscle activation patterns.5,19 This leads to the first objective of this study to use an EMG-driven model to estimate ATT, anterior shear forces, and knee ligament loading during ACL-deficient gait.

Anterior–posterior knee stability has been shown to be related to posterior inclination angle of the tibial plateau in previous experiments. An in vivo radiographic study showed increase in ATT during weight bearing for both healthy and ACL-deficient knees.8 Two in vitro studies found that the resting position of the tibia shifted anteriorly as posterior tibial slope increased in healthy knees, and the in situ ACL force and ACL strain decreased.11,12 The limitation of these in vitro studies is that the loads added to the cadaveric knees were smaller than those of in vivo weight-bearing conditions during gait. This leads to the second objective of this study to use a numerical model to study the influence of increasing posterior tibial slope on ATT and ligament loading.

In this article, we will describe for the first time an EMG-driven model that incorporates a knee–ligament model, and apply this approach to estimate ATT, anterior shear forces, and ligament loading in the knee joint of a healthy subject and an ACL-deficient subject. We validated this model by comparing with previous published results, and used the calculated results to discuss how the ACL-deficient subject compensated for the loss of the ACL. We will also compare the difference of ATT and ligament loading between two knees that have different tibial slope inclinations, and study the influence of increasing tibial slope on knee biomechanics.

Methods

Data Collection

One male healthy subject (mass 60.5 kg, height 1.70 m) and one male patient without an ACL (mass 74.0 kg, height 1.72 m; right leg was the affected leg) gave informed consent before participating in this study. The experimental protocol was approved by the Human Subjects Review Board of the University of Delaware. EMG, joint position, and force plate data were collected from four walking trials in our motion analysis lab. In addition, four isometric maximal voluntary trials and four isokinetic maximal voluntary trials (60°/s) of knee flexion/extension on a Biodex dynamometer (Shirley, NY, USA), five straight running trials (4.0 m/s), and five running sidestepping trials (to 45° from the direction of travel, 4.0 m/s) were collected. EMGs were collected from nine muscles of the right knee using surface electrodes, including rectus femoris (RF), VL, vastus medialis (VM), semimembranosus (SM), biceps femoris long head (BFL), medial gastrocnemius (MG), and LG. In this study we simulated the stance phase of the right leg during the walking trials.

Computational Model

We calculated ATT and ligament loading through a two-step procedure.

EMG-Driven Model

First, an EMG-driven model4,29,38 was used to estimate the muscle forces of the right leg. EMG and kinematic data were used as inputs. The model included the following four main parts.

Anatomical Model

The lower limb anatomical model was developed using SIMM (Musculographics Inc.),10 and scaled based on the subject’s static marker measurements using Opensim.9 The model included ten musculotendon actuators crossing the knee joint. These were RF, VL, VM, vastus intermedius (VI), SM, semitendinosus (ST), BFL, biceps femoris short head (BFS), MG, and LG.

Muscle Activation Dynamics Model

The muscle activation dynamics model transformed raw EMG to muscle activation. Firstly, the raw EMG were high-pass filtered, then full wave rectified, normalized by peak rectified EMG values obtained from maximum voluntary contraction trials, and then low-pass filtered.4 In the second step, we used a discrete form of a linear second-order differential equation to process the filtered EMG, and took the electromechanical delay into account.44 Finally, reported nonlinear EMG–force relationships were accounted for by using a logarithmic function for low levels of EMG and a linear function for high values to transform the processed EMG to the muscle activation.30 For the three muscles that did not have surface EMG measurements, the muscle activations were calculated using the same procedure as Lloyd and Besier.29

Muscle Contraction Dynamics Model

The muscle contraction dynamics model employed a modified Hill-type muscle model to calculate individual muscle forces. The muscle activations and joint kinematics were used as inputs. The muscle–tendon unit was modeled as a muscle fiber in series with a tendon. The muscle fiber had a contractile element13,49 in parallel with an elastic element36 and a damping element.37 Once muscle forces were computed, they were multiplied by respective muscle moment arms and summed to determine total joint moments.

Calibration Process

Since some of the model parameters (such as muscle activation model parameters, optimal fiber length, tendon slack length, electromechanical delay, and strength factor of muscles) were difficult to get from in vivo measurements, we used a calibration process to obtain a set of model parameters for each subject that could accurately estimate the knee joint moments during a variety of dynamic tasks. It was assumed that the calculated EMG-driven joint moments should equal the inverse dynamic calculated joint moments. An optimization algorithm20 was employed, and the parameters in the EMG-driven model were tuned to minimize the difference between the EMG-driven knee joint moments and inverse dynamics knee joint moments. The calibration trials included five trials: (1) maximal isokinetic knee flexion, (2) maximal isokinetic knee extension, (3) walking, (4) straight running, and (5) running sidestepping. All the model parameters were constrained within a reasonable range as suggested in the literature.29

After the parameters were calibrated, the tuned EMG-driven model was used to predict knee joint moments for the 17 other trials for each subject. Since direct measurement of in vivo muscle force was difficult and invasive, it is impractical to validate numerical models directly using muscle force measurements. Therefore, joint moments measured using an inverse dynamics approach are typically used as an indirect way to validate EMG-driven models.29 In this study, we validated the tuned EMG-driven model by comparing the predicted joint moments with inverse dynamics joint moments. The coefficient of determination (R 2) between the predicted joint moments and the inverse dynamics joint moments, the normalized root mean square (RMS) error (normalized to peak-to-peak joint moment), and the normalized maximal error were calculated for each prediction trial. Once the tuned EMG-driven model had been validated, the model was then used to estimate the muscle forces during the walking trials, and the calculated knee muscle forces were used as inputs to the following knee joint model to calculate ATT and ligament forces.

Knee Joint Model

A knee model that incorporated a knee–ligament model was developed in the second step. The model included the tibiofemoral joint and the patellofemoral joint in the sagittal plane (Fig. 1a), and three segments: the femur, tibia, and patella. The origin of the femur coordinate system was located at the midpoint between medial and lateral femoral epicondyles, with the x-axis pointing posteriorly and the y-axis pointing proximally along the long axis of the femur. The origin of the tibia coordinate system was located at the center of mass of the tibia, with the u-axis pointing posteriorly and the v-axis pointing proximally along the long axis of the tibia.

FIGURE 1
figure 1

(a) The two-dimensional patellofemoral and tibiofemoral joint models. Knee joint angle (θ) represented the flexion/extension of the tibia coordinate system relative to the femur coordinate system. Muscle forces (F VAS and F RF), patellar tendon force (F PT), and patellofemoral contact force (F pf) were applied on the patella. F pf was applied along the normal direction of internal edge of the patella at patellofemoral contact point, P c. α was the angle between the long axis of the patella and the y-axis; δ VAS, δ RF, and β were the angles between the muscle or ligament lines and the y-axis. (b) The forces applied on the tibia included: hamstrings, F hs (including SM, ST, BFL, and BFS), and gastrocnemius, F gas (including MG and LG), calculated from the EMG-driven model in the first step, external forces of the ankle joint (F ax and F ay) calculated from the inverse dynamics in the first step, patellar tendon force (F PT) calculated from the patellofemoral joint model, ligament forces (F ACL, F PCL, F MCL, and the forces of the other ligaments, F others) from the ligament model, inertia forces, gravity, and tibiofemoral joint contact force (F contact)

Full size image
The Patellofemoral Joint Model

The approach used to model the patellofemoral joint in this study was similar to that used by Shelburne and Pandy,39 and Liu and Maitland27 (Fig. 1a). The patella was approximated as a rectangle,47 and the length and width data were taken from magnetic resonance imaging measurement.48 The femoral contour of the patellofemoral joint was adapted as an ellipse.27 The geometrical parameters of the patellofemoral joint were summarized in Table 1. The origin and insertion points of muscles and patellar tendon were from Delp et al. 10 All the geometrical parameters in the knee joint, including insertion points of muscles and ligaments, segment lengths, and femoral contour sizes in the knee joint model, were scaled based on the subject’s static marker measurements using Opensim.9 Muscle forces (RF, VL, VM, VI) calculated from the first step were applied on the patella along the muscle lines. The contact between the femoral condyles and the internal edge of the patella was assumed to be rigid and frictionless. We assumed that the mass of the patella was negligible in comparison to other segments, and the patellar tendon was inextensible at all joint angles. Under these assumptions, the three force and moment equilibrium equations of the patella could be combined into just one equation.39 The orientation of the patella, α, was the only unknown and could be solved iteratively. Once α was known, the patellar tendon force, F PT, could be solved through the force equilibrium equations.

Table 1 Parameters used in the knee model
Full size table
The Tibiofemoral Joint Model

The tibiofemoral joint had three degrees-of-freedom, including knee joint angle (θ), anterior–posterior position (x t), and superior–inferior position (y t) of the tibia coordinate system’s origin relative to the femur coordinate system (Fig. 1a). The contour of the tibial plateau was described as a straight line. Since posterior inclination angle of the tibial plateau was reported to be highly variable between populations,32 we chose 4°16 and 8°31 as the two tibial slopes relative to the u-axis in this study. The femoral contour of the tibiofemoral joint was adapted as an ellipse.47 The parameters of the two contours were described in Table 1. The contact between the femoral condyles and tibial plateau was assumed to be rigid and frictionless. Under these assumptions, y t could be represented by a function of x t and θ.1 Since knee joint angle was determined from measurement, x t was the only unknown degree-of-freedom under the geometric constraints.

The Knee–Ligament Model

A knee–ligament model was developed to calculate ligament forces. The cruciate and collateral ligaments, the posterior capsule, and the anterolateral structures of the knee were modeled using 14 elastic bundles, and each bundle was represented by a straightline from its origin to its insertion point.33,40 They were the anterior and posterior bundles of the ACL, the anterior and posterior bundles of the posterior cruciate ligament (PCL), the anterior, central, and posterior bundles of the superficial MCL, the anterior and posterior bundles of the deep MCL, the lateral collateral ligament, the popliteofibular ligament, the medial and lateral posterior capsule, and the anterolateral structures. Each ligament bundle was modeled as a nonlinear elastic element, described by a nonlinear force–strain relationship.3 For the ACL-deficient knee, the ACL force was set to be zero. The model also included a posterior constraint force contributed by the medial meniscus, which increased linearly with respect to the ATT.40

The stiffness and zero load length of the ligament bundles, and the point at which the constraint force of the medial meniscus begins and the increasing rate were tuned so that the calculated ATTs of the knee–ligament model in response to a 100-N anterior force matched the ATT results of previous in vitro experiments on healthy knees and ACL-deficient knees.25,26,43 A program coded in Matlab (The MathWorks, Inc., Natick, MA, USA) was developed, and lsqnonlin() with ‘trust-region-reflective’ as the optimization algorithm was chosen to solve for the ligament parameters. All the ligament parameters were given an initial value as suggested by the literature,33,40 and they were allowed to change within ±20% to account for individual variation. First the knee–ligament model was used to calculate the resulting ATTs by adding a 100-N anterior force, and then the knee–ligament model without the medial meniscus (set the medial meniscus force to be zero) was used to calculate the resulting ATTs. The differences between the ATTs of these two calculations (with/without the medial meniscus) and the corresponding experimental ATT data (with/without the medial meniscus) were summed as the cost function. We also conducted a sensitivity study on this tuning process. Another two sets of ATT data were used as inputs to match the calculated results with during the optimization. The first one used the experimental data plus 3 mm, and the second one used the experimental data minus 3 mm. The percentage change of the ligament parameters and the changes in the resulting ATT and ligament loading were calculated to find out how our model responded to different input ATTs.

The Tibia Equilibrium Model

All forces applied on the tibia (Fig. 1b) were assumed to be balanced at each time step during stance phase, and the foot was not included in this model. The inertia forces caused by rectilinear acceleration of the tibia were calculated using the measured rectilinear acceleration of the center of mass, and the inertia forces caused by rotation of the tibia, including Coriolis force, centrifugal force, and Euler force, were ignored since they were small during the stance phase of gait.15 Gravity of the tibia was applied on the center of mass of the tibia. The tibiofemoral joint contact force (F contact) was applied on the tibiofemoral contact point along the normal direction of the tibial plateau. The magnitude of tibiofemoral joint contact force was unknown. Anterior–posterior position of the tibia (x t) was another unknown, and it determined the relative positions of the femur to the tibia (x t, y t, θ), which were used as inputs to the patellofemoral joint and tibiofemoral joint model. At each time step two force equilibrium equations of the tibia in u and v directions were obtained. Two independent unknown variables, including the position variable (x t) and the magnitude of tibiofemoral joint contact force (F contact), were solved through iterations until the force equilibrium of the tibia was satisfied. A program coded in Matlab (The MathWorks, Inc.) was used to solve the nonlinear equations, and fsolve() with ‘trust-region-dogleg’ as the optimization algorithm was chosen as the solver. After x t and F contact were calculated, the ATT, anterior shear forces, and ligament loading could be calculated using the patellofemoral joint model, tibiofemoral joint model, and knee–ligament model.

Results

Validation of the Tuned EMG-Driven Model

Table 2 showed the prediction results using the tuned EMG-driven model for both subjects. The average R 2 values were 0.848 and 0.826, and the normalized RMS errors and maximal errors were small, so that the predicted EMG-driven joint moments were close to the inverse dynamics calculated joint moments in a variety of tasks. This validated our tuned EMG-driven model.

Table 2 Average coefficient of determination (R 2) between the predicted EMG-driven joint moments and the inverse dynamics calculated joint moments, the normalized RMS error (normalized to peak-to-peak joint moment), and the normalized maximal error during 17 prediction trials
Full size table

Validation of the Knee Model

The calculated ATT patterns of the healthy subject (Fig. 3a) was similar to previous in vivo experiments using six degrees-of-freedom goniometer24,50 and dual fluoroscopic imaging technique23 (Fig. 3b). The calculated ATT patterns in the ACL-deficient knee (Fig. 3a) was also consistent with previous in vivo experiments24,50 (Fig. 3b). The increase in ATT after ACL rupture has also been shown in a number of in vitro 2,14,22,25,26,42,43,46 studies (Fig. 3b). The average RMS error between the calculated ATTs and the experimental ATTs was 3.22 mm for the healthy knee with 4° tibial slope and 3.23 mm for the ACL-deficient knee with 4° tibial slope, and the differences were within the standard deviation of the measurements.24

Our results of the healthy knee with 4° tibial slope showed that the peak anterior shear force applied by tibiofemoral joint contact force was 0.30 BW, and Shelburne et al. 40 generated result of 0.11 BW using a numerical model. Our calculation was close to two in vivo studies that measured the peak anterior shear forces (0.2918 and 0.30 BW7) using an instrumented knee implant during walking. The main passive restraint of anterior shear force was the MCL in the ACL-deficient knee as predicted by our model (Figs. 4a and 4b). Previous in vitro experiments showed similar results that the in situ force of the MCL would increase after transection of the ACL.22 Meniscus also contributed to the balance of anterior shear force, and this is consistent with previous in vitro experimental result.2

The magnitudes of ligament force estimation in our model were comparable to previous numerical results. The maximum ACL force in the healthy knee of our model (0.37 BW for the knee with 4° tibial slope and 0.54 BW for the knee with 8° tibial slope) was similar to the previous numerical models15,40 (Table 4). The maximum MCL force in the ACL-deficient knee was three times bigger than that of the healthy knee in our model, and this is consistent with the numerical result of Shelburne et al. 40 (Table 4). Our calculated maximum MCL force in the ACL-deficient knee (0.59 BW for the knee with 4° tibial slope and 0.88 BW for the knee with 8° tibial slope) was bigger than the numerical result of Shelburne et al. 40 (0.17 BW).

Ligament Parameter Adjustment

After adjustment of the ligament model parameters, the calculated anterior–posterior laxity of the healthy knee and the ACL-deficient knee (with/without medial meniscus) were within the range of in vitro measurements (Fig. 2). It gave us an adjusted ligament model that could predict the response of the knee ligaments to external loads. Table 3 showed that by increasing or decreasing the in vitro ATT data by 3 mm the percentage change of the tuned ligament parameters was 0–8%. The RMS errors of the ATT change were relatively small (0.48–2.36 mm), and specially for the ACL-deficient subject the R 2 values of the ATT change were greater than 0.96, which indicated a small change in the general pattern of the ATTs. The RMS errors of both the ACL force and MCL force were 54.2–212.8 N, and the R 2 values were 0.75–0.97, which indicated a small change in the general force patterns.

FIGURE 2
figure 2

Calculated ATT for (a) the healthy knee, (b) the ACL-deficient knee, and (c) the ACL-deficient knee without meniscus in response to a 100-N anterior force using the adjusted ligament model (black line). Also shown for comparison are in vitro experimental data reported by Levy et al.,25,26 Haimes et al.,14 and Sullivan et al. 43 The vertical lines on the experimental data represent plus and minus one standard deviation. We chose 0° to 45° as the data of interest because the knee joint angles of these two subjects during stance phase were within this range

Full size image
Table 3 The differences between the new results using the in vitro ATT data minus 3 mm and the in vitro ATT data plus 3 mm as inputs and the results using the original in vitro data as inputs during tuning of the knee-ligament parameters
Full size table

Anterior Tibial Translation

Here, we report the results of one walking trial of the healthy knee and the ACL-deficient knee using our model. The ACL-deficient knee had increased ATT compared to the healthy knee throughout stance phase for both knees with 4° and 8° tibial slope (Fig. 3). The ATT profiles of both healthy and ACL-deficient knees peaked twice at about 20% of stance phase and near toe off. The ATT of the knee with 8° tibial slope was greater than that of the knee with 4° tibial slope for both knees, which showed increasing tibial slope in our model would increase ATT.

FIGURE 3
figure 3

(a) Calculated ATT in the healthy and ACL-deficient knee with 4° and 8° tibial slopes during stance phase using our model. The ATT was measured as the anterior translation of the tibia relative to the midline between medial and lateral femoral epicondyles. (b) Adapted from previous in vivo measurements of ATT23,24,50

Full size image

Anterior Shear Force

In the healthy knee, the ACL was the major passive restraint to balance anterior shear force (Fig. 4a). Tibiofemoral joint contact force and patellar tendon force were the major sources of anterior shear force in the healthy knee, and they peaked at about 20% of stance phase. The anterior shear force profiles of the healthy knee with 8° tibial slope were similar to Fig. 4a, except that the anterior shear force generated by tibiofemoral joint contact force and the main passive restraint force generated by the ACL were bigger for 8° compared with those of 4° (not shown in the figure). The peak anterior shear force generated by tibiofemoral joint contact force was 0.28 BW for 4° tibial slope and 0.58 BW for 8° tibial slope, and the peak constraint force generated by the ACL was 0.22 BW for 4° tibial slope and 0.35 BW for 8° tibial slope.

FIGURE 4
figure 4

Anterior shear forces applied to the tibia by muscles, ligaments, external force, and joint contact force in (a) the healthy knee with 4° tibial slope and (b) the ACL-deficient knee with 4° tibial slope. The ACL provided more than 90% of the passive restraint provided by soft tissues in the healthy knee, and the MCL provided more than 90% of the passive restraint provided by soft tissues in the ACL-deficient knee, so the other ligaments were not shown. BW, body weight

Full size image

In the ACL-deficient knee, the MCL was the major passive restraint of anterior shear force instead of the ruptured ACL (Fig. 4b). Meniscus and posterior capsule also contributed to balance anterior shear force (data not shown in the figure). We found strong co-contractions of hamstrings during early stance phase and strong co-contractions of gastrocnemius throughout stance phase. Both hamstrings and gastrocnemius functioned as active restraints of anterior shear force and peaked during early stance phase. Tibiofemoral joint contact force and patellar tendon force were the major sources of anterior shear force in the ACL-deficient knee, and they peaked at about 20% of stance phase. The anterior shear force profiles of the ACL-deficient knee with 8° tibial slope were similar to Fig. 4b, except that the anterior shear force generated by tibiofemoral joint contact force and the main passive restraint force generated by the MCL were bigger for 8° compared with those of 4° (data not shown in the figure). The peak anterior shear force generated by tibiofemoral joint contact force was 0.32 BW for 4° tibial slope and 0.67 BW for 8° tibial slope, and the peak constraint force generated by the MCL was 0.32 BW for 4° tibial slope and 0.40 BW for 8° tibial slope.

Ligament Loading

In the healthy knee, the ACL was loaded for most of the stance phase, and the MCL, lateral collateral ligament, and posterior capsule were the other three loaded ligaments (Figs. 5a and 5b). The ACL force had two peaks during early stance phase and mid-stance phase. The ligament forces of the healthy knee with 8° tibial slope were greater compared to those of the healthy knee with 4° tibial slope (Table 4). In the ACL-deficient knee, the MCL was loaded for most of the stance phase, and the meniscus and posterior capsule were also loaded at a small magnitude (Figs. 5c and 5d). The lateral collateral ligament and PCL were not loaded at all during stance phase. The ligament forces of the ACL-deficient knee with 8° tibial slope were greater compared to those of the ACL-deficient knee with 4° tibial slope (Table 4).

FIGURE 5
figure 5

Ligament loading during stance phase in (a) the healthy knee with 4° tibial slope, (b) the healthy knee with 8° tibial slope, (c) the ACL-deficient knee with 4° tibial slope, and (d) the ACL-deficient knee with 8° tibial slope. The unloaded ligaments were not shown here

Full size image
Table 4 Comparison of the maximal ligament forces during stance phase between the current study and two previous simulation studies
Full size table

Discussion

We have developed a biomechanical model that used EMGs as input and incorporated a knee–ligament model to estimate ATT and ligament loading during healthy and ACL-deficient gait. The calculated results were similar to previous experimental and numerical studies, and provided insights on how the ACL-deficient subject adapted their gait following the loss of the ACL.

A unique feature of EMG-driven model is that it could account for the abnormal EMG patterns and muscle activation strategies that the subject without an ACL was using, and it helped reveal the subject’s muscle compensation strategies for the loss of ACL. This is what makes our approach different from the previous numerical models,27,40 and why the approach is more suitable for the study of pathological gait. An in vivo experimental study showed that the gastrocnemius muscles of the ACL-deficient knees were active during stance phase, and hamstrings were more active during early stance, but the ATT increased more rapidly in spite of the strong co-contraction.24 In our study, we also found strong co-contraction of gastrocnemius throughout stance phase and longer duration of hamstrings activity during early stance phase in the ACL-deficient knee. The calculated results showed that the tibia of the ACL-deficient knee shifted more anteriorly throughout the stance phase. This suggested that the ACL was the main restraint of anterior shear force in the knee, and it could not be compensated using knee flexors in some of the ACL-deficient subjects.

We validated our model by comparing previous experimental and numerical studies, and this gave us confidence on the accuracy of the model calculations. The calculated ATT patterns of both subjects were similar to previous in vivo experiments23,24,50; the calculated anterior shear force applied by joint contact force in the healthy knee was close to previous in vivo experiments7,18; the calculated ligament forces were comparable to previous numerical models.15,40 Since it is difficult to finish in vivo muscle force or ligament force measurements, we did not validate our calculated muscle forces and ligament forces by directly comparing to measurements. In the future, we may use in vivo imaging measurements as an alternative. We can use the measurements of ligament bundle length as references to validate our ligament model, and use the measured muscle fiber lengths to validate our EMG-driven muscle model.

The calculated results could be used to explain knee biomechanics during healthy and ACL-deficient gait. We found two ATT peaks in both the healthy and ACL-deficient knees occurring during early stance phase and near toe off. The two ATT peaks coincided with the two knee extension torque peaks. A bigger patellar tendon force during peak of extension torque would add more anterior shear load to the tibia, and this would lead to an increased ATT. The main passive restraint of anterior shear force was the ACL in the healthy knee and the MCL in the ACL-deficient knee. Hamstrings and gastrocnemius functioned as active restraints of anterior shear force in the ACL-deficient knee. The knee flexor forces peaked during early stance phase, which coincided with the maximum ATT and anterior shear force.

Tibiofemoral joint contact force and patellar tendon force were found as the main sources of anterior shear force for both healthy and ACL-deficient knees, and the tibial slope is a crucial parameter in determining knee mechanics. Our calculated anterior shear force applied by tibiofemoral joint contact force in the healthy knee was bigger (~0.30 BW for 4° tibial slope; ~0.60 BW for 8° tibial slope) than the result of Shelburne et al. 40 (~0.11 BW). They used 2° and 7° as the posterior inclination angles of the lateral and medial tibial plateaus, respectively, whereas we used the same angle for both medial and lateral sides in our knee model, with 4° and 8° tibial slopes. Since tibiofemoral joint contact force is along the normal direction of the tibia plateau, the anterior component of tibiofemoral joint contact force is determined by the sine of the posterior inclination angle (sin 2° and sin 7° vs. sin 4° and sin 8°), which led to the big difference between the two results. In literature the posterior inclination angle of the tibial plateau has been reported to have a mean value between 4° and 14°.32 Two in vivo studies showed that the peak anterior shear forces during walking were 0.2918 and 0.30 BW,7 and it was close to our results of 4° tibial slope. Therefore, we think our model was more valid, and using a subject-specific measurement would be a more accurate way in future studies. The bigger anterior shear force generated by tibiofemoral joint contact force in our model could also explain the higher estimation of the MCL force in our model compared to that of Shelburne et al.,40 because a bigger MCL force was required to balance the bigger anterior shear force.

The results of our two knee models with 4° and 8° posterior tibial slope showed that increasing the tibial slope in our model would increase the ATT and ligament loading for both healthy and ACL-deficient knees. This has also been shown in an in vivo radiographic study8 and two previous numerical studies.21,28 During weight-bearing conditions tibiofemoral joint contact force generates an anteriorly directed force component, and this anterior component increases as the tibial slope increases. This increased anterior component can produce a corresponding anteriorly directed translation of the tibia and increase ligament loading that is necessary to balance the anterior shear force. The positive relationship between posterior tibial slope and ACL injury rate has been verified in a case–control study.17 Females in general have a steeper tibial slope compared to males,16 and the bigger anterior shear force generated by tibiofemoral joint contact force during weight-bearing conditions moves their tibia more anteriorly and exposes them to an increased risk of ACL injury compared to males.

The sensitivity study of the knee–ligament model showed that the calculated results would be influenced by the inputs of in vitro ATT data during the tuning of the ligament parameters. The changes in the calculated ATTs were relatively small, especially for the ACL-deficient subject. The ligament force patterns were not changed much after changing the in vitro ATT inputs to the ligament tuning model, but the force magnitudes were changed by an amount of 100 or 200 N. This study showed that it is important to have subject-specific in vivo ATT data in order to tune the knee–ligament parameters more accurately.

In future studies we will incorporate imaging measurements to make our model more subject-specific, and use the model to provide useful information for patient-specific surgery procedure and rehabilitation programs. Kvist and Gillquist24 found the inter-subject variance of maximal ATT during walking to be ±2.3 mm for healthy knees and ±4 mm for ACL-deficient knees, and the ligament loading would also vary across different subjects. This big variability between subjects would require a more subject-specific model to identify the individual difference in our future model development. First, we will include imaging measurements to provide subject-specific anthropometric parameters, including the parameters of the femur and tibia contours, the muscle and ligament attachment points. Second, we will engage in vivo magnetic resonance imaging41 or ultrasound6 measurements during muscle contractions to make our EMG-driven muscle model more subject-specific. Third, we will use in vivo ATT data from MRI or X-ray measurements in response to different external loads to tune our ligament model’s parameters, and this will make our ligament model more subject-specific. The increase in the model accuracy can strengthen our model as a tool to help decisions on the surgery and rehabilitation. For example, we can use our model to predict the ligament graft loading after different surgery procedures (patellar tendon graft vs. hamstring tendon draft, graft tunnel position, initial graft tension, etc.), and find the best procedure that would restore the injured knee to normal function. We can also use our model to calculate the ligament loading during various activities of rehabilitation, such as strengthening program and running program, and we will examine if the rehabilitation protocol is safe for a specific patient without reinjuring the ligaments.

There are some limitations in this study. First, this model was confined to the sagittal plane, and in future studies we will expand our model to the frontal and transverse plane. Second, we assumed that the joint contact was rigid, and ignored the deformability of the meniscus and cartilage. Our calculated ATT results of both healthy and ACL-deficient knees were close to in vivo measurements,23,24,50 and an in vitro study26 has shown that the anterior–posterior laxity of healthy knees would not be affected after removal of the meniscus, so we think that the ATT results in our study were not affected by the rigidity assumption. The ligament force estimations were likely to be higher because a deformable joint contact was not taken into account.39 Third, experimental data were taken from a minimal amount of subjects in this study. In the future, this approach will be applied to a patient population for statistical comparison.

This is the first time that an EMG-driven model has been incorporated with a knee–ligament model to estimate ligament loading in a healthy knee and an ACL-deficient knee. The calculated results were consistent with previous experimental and numerical studies, and this gave us confidence that our model could be used to study how ACL-deficient knees compensate for the loss of the ACL. We found increase in ATT throughout stance phase for the ACL-deficient knee compared to the healthy knee. The MCL was the main passive restraint to anterior shear force in the ACL-deficient knee, and knee flexors were used as active restraint to help balance anterior shear force. We anticipated that increasing the tibial slope would increase the resulting ATT and ligament forces in both healthy and ACL-deficient knees based on the results of our model.