Simulation of anterior cruciate ligament deficiency in a musculoskeletal model with anatomical knees.

Abnormal knee kinematics and meniscus injury resulting from anterior cruciate ligament (ACL) deficiency are often implicated in joint degeneration even though changes in tibio-femoral contact location after injury are small, typically only a few millimeters. Ligament reconstruction surgery does not significantly reduce the incidence of early onset osteoarthritis. Increased knowledge of knee contact mechanics would increase our understanding of the effects of ACL injury and help guide ACL reconstruction methods. Presented here is a cadaver specific computational knee model combined with a body-level musculoskeletal model from a subject of similar height and weight as the cadaver donor. The knee model was developed in the multi-body framework and includes representation of the menisci. Experimental body-level measurements provided input to the musculoskeletal model. The location of tibio-menisco-femoral contact as well as contact pressures were compared for models with an intact ACL, partial ACL transection (posterolateral bundle transection), and full ACL transection during a muscle driven forward dynamics simulation of a dual limb squat. During the squat, small changes in femur motion relative to the tibia for both partial and full ACL transection push the lateral meniscus in the posterior direction at extension. The central-anterior region of the lateral meniscus then becomes "wedged" between the tibia and femur during knee flexion. This "wedging" effect does not occur for the intact knee. Peak contact pressure and contact locations are similar for the partial tear and complete ACL transection during the deep flexion portion of the squat, particularly on the lateral side. The tibio-femoral contact location on the tibia plateau shifts slightly to the posterior and lateral direction with ACL transection.

INTRODUCTION

It is estimated that 27 million adults in the United States have clinical osteoarthritis with knee osteoarthritis being the most prevalent form affecting 28% of adults over age 45 and 37% of adults over age 60 [1]. Although the prevalence and debilitating nature of osteoarthritis is well documented, the etiology of this chronic disease is not completely understood. The most significant cause of osteoarthritis in young adults is joint trauma including tears of the menisci or ligaments. For example, 12 years after rupture of the anterior cruciate ligament (ACL), 51% of a female population with a mean age of 31 years had radiographic evidence of osteoarthritis [2]. Similarly, 43% of patients with intact ACLs and isolated limited meniscectomy had radiographic evidence of osteoarthritis after 16 years [3]. In addition, the severity of osteoarthritis correlates with the severity of meniscal injury and the amount of tissue removed [4].

Abnormal knee kinematics and meniscus injury resulting from ACL deficiency are often implicated in joint degeneration [5] while ACL reconstruction surgeries attempt to restore normal tibio-femoral motion and prevent long-term consequences [6]. The observed differences in tibio-femoral kinematics between normal and ACL deficient knees are small, typically only a few millimeters [7, 8]. Studies have also shown that ACL reconstruction does not significantly reduce the incidence of early onset osteoarthritis [2]. Increased knowledge of tibiofemoral joint contact mechanics would increase our understanding of the effects of ACL injury and help guide ACL reconstruction methods. Medical imaging has been used to quantify motion and cartilage contact in normal and ACL deficient knees including knee loading in a magnetic resonance imaging (MRI) tunnel [9] and MRI combined with fluoroscopic images captured during a weight-bearing lunge [10]. Computational models capable of predicting tibio-menisco-femoral contact mechanics during movement would also provide a valuable tool for increasing our understanding of normal joint mechanics as well as knee injury and joint degeneration.

The main objective of this study was to develop a muscle driven forward dynamics model of the lower extremities in a dual limb squat, validate the computational model against experimental data and then simulate the same motion in two different ACL deficient conditions. The model combined a cadaver anatomical knee with anthropometric and motion data from a female subject of similar height and weight as the cadaver donor. The model was developed in the multibody framework and it was validated against measurements of the ground reaction forces and surface electromyography (EMG) taken during the squat motion. The multibody model included representation of the menisci with an intact ACL, transected posterolateral ACL bundle, and full ACL transection. The contact pressures on the tibia plateau for the intact and ACL deficient cases were compared during a muscle driven forward dynamics simulation.

MATERIALS AND METHODS

Computational Knee Model

A previously developed anatomical computational knee model was used for this study. The multibody model was created in MD Adams (MSC Software Corporation, Santa Ana, CA) and has been previously described [11]. Only a brief summary is provided here. Knee geometries (tibia, femur, articular cartilage and menisci) were derived from magnetic resonance images (MRI) of a fresh frozen cadaver knee (55 years old, female, left knee, 170 cm height, 72 kg mass). After MRI, the cadaver knee was mounted in a dynamic knee simulator (Kansas Knee Simulator, University of Kansas, Lawrence, KS) and manipulated through a walk cycle. During testing, the forces applied by the actuators of the machine were recorded and bone motion was measured using rigid body markers attached to the femur, tibia, and patella with an Optotrak 3020 system (Northern Digital Inc., Waterloo, Ontario). For calculation of ligament zero-load lengths, the motion of the tibia and femur rigid body markers was recorded while the joint was manually moved through its full range of motion with a minimum force applied (as judged by the experimenter). After testing, the knee was disarticulated and the insertion sites for the anterior cruciate ligament (ACL), posterior cruciate ligament (PCL), medial collateral ligament (MCL), lateral collateral ligament (LCL), patellar tendon, menisci horn attachments, and menisci transverse ligament were measured with an Optotrak digitizing probe.

In the computational model, ligament bundles were represented as one-dimensional non-linear springs. The model included two bundles for the ACL [12] and PCL [13] and three bundles for the MCL and LCL [14]. Non-linear splines were used to describe the force-displacement curve of each ligament including the non-linear “toe” region. The splines were derived from the ligament force as a function of strain, the length of each ligament in the position it was constructed, and the measured zero-load length. The zero-load length of each bundle was determined by calculating the maximum straight-line distance between insertion and origin sites throughout the experimentally measured full range of motion and then applying a correction percentage [11]. The force-length relationship for each ligament is described by: [15, 16]

where k is a stiffness parameter, ε is a spring parameter assumed to be 0.03, ε is the engineering strain of each ligament bundle, l is ligament bundle length, and l is the zero-load length. Values of k for each ligament bundle came from Wismans et al. and Blankevoort et al. [15, 16]. Each 1-D spring also included a parallel damper and a damping coefficient of 0.5 Ns/mm was used for each ligament bundle.

As described in Guess, Liu et al., [11], the medial and lateral tibia plateau cartilage geometries were divided into multiple hexahedral rigid bodies. The size of each hexahedral element’s cross-section in the transverse plane of the tibia was 4 x 4 mm and a total of 61 lateral and 72 medial elements were created. Each cartilage element was connected to tibia bone with a fixed joint located at the center of each tibia cartilage-bone interface. A deformable contact constraint was defined between each tibia cartilage element and the femur cartilage geometry. The contact model used for all articulating surfaces in the knee was defined as:

where F is the contact force, δ is the interpenetration of geometries, is the velocity of interpenetration, k is a spring constant, n is the compliance exponent, and B(δ) is a damping coefficient. The contact parameters defined between each cartilage rigid body and the femur cartilage were derived by optimizing the contact parameters to match contact pressure predictions of a finite element model [11]. Specifically, the contact parameters were optimized to minimize the difference between pressure predictions of the multibody model and that of an identically loaded and constrained finite element model. The resulting parameters were k = 140 N/mm1/1.3, n = 1.3, and B = 5 Ns/mm.

Multibody models of the menisci were created by radially sectioning the lateral and medial menisci geometries into 17 elements each and assigning mass properties to every segment based on its volume and a density of 1100 kg/m3 [17]. The meniscus rigid body elements were connected to neighboring elements by the following stiffness matrix:

where Fθ,r,z and Tθ,r,z are the translational and torsional forces between elements acting in the circumferential, radial, and axial directions. Kθ, Kr, Kz, Kθr, Kθz, Krz, Tθ, Tr, and Tz are the stiffness matrix parameters, θ, r, and z are relative translational displacements and a, b, and c are relative rotational displacements. The stiffness matrix parameters are shown in Table . Values for the stiffness parameters were derived through an optimization process that minimized the displacement error between identically loaded finite element and multibody menisci models [18].

Each meniscus was attached to the tibia via 4 horn attachments (2 posterior and 2 medial per meniscus) modeled as simple tension only springs with a spring constant of 1000 N/mm [19]. The transverse ligament was modeled as a tension only spring with a stiffness of 200 N/mm. A parallel damper with damping coefficient of 0.5 Ns/mm was included for each horn attachment and the transverse ligament spring. Deformable contacts using Eq. 3 were defined between each meniscus element and the femur cartilage. In addition, a deformable contact (k = 80 N/mm1/2.5, n = 2.5) was defined between each meniscus element and every tibia cartilage element that might come in contact with the meniscus element during movement (122 contacts medial side and 134 contacts lateral side). A deformable contact was also defined between the patella cartilage and femur cartilage geometries (Fig. ). To replicate the anatomical knee force-displacement relationships for the right leg, a right knee was created by mirroring the geometries and ligament insertion sites of the left knee. As the focus of the study was on the cadaver based left knee, the menisci were not included in the mirrored right knee.

Gait Measurements

With informed consent, motion, ground reaction forces, and surface electromyography (EMG) were collected on a young female subject (age 27) of similar height (175 cm) and weight (64 kg) as that of the cadaver knee donor. Markers were placed on the subject according to the plug-in-gait protocol and a dual-limb squat activity was performed where the subject’s right and left feet were isolated on different force plates. During the squat the heels were always in contact with the ground and the subject was instructed to start with the knee slightly flexed. A squat activity with an initial flexion angle matching the flexion of the computational model was chosen (Fig. ). The squat began with a knee flexion angle of 16° (0.25 s, extension) and the deepest flexion angle was 87° (1.625 s, flexion).

Musculoskeletal Model

LifeMOD (LifeModeler, Inc., San Clemente, CA) was used for development and simulation of the musculoskeletal model (Fig. ). LifeMOD is a virtual human modeling and simulation software add-on to MD ADAMS. The program is used for development of musculoskeletal structures and dynamic movement simulations based on subject anthropometrics and captured motion data. For forward dynamics simulations, LifeMOD predicts the muscle forces that will reproduce measured body motion. The model developed for this project included the pelvis, right lower limb, and left lower limb. The right and left hip were modeled as spherical joints as well as the right and left ankle. Anthropometric measurements from the test subject were used to scale the model. The previously developed left knee and mirrored right knee were placed in the musculoskeletal model with the femur and tibia of the anatomical knees rigidly attached to the musculoskeletal model femur and tibia. Measured motion data from the gait testing was imported into the model and knee flexion-extension motion was used to position the anatomical knees relative to the musculoskeletal model. This was accomplished by systematically repositioning the anatomical knees until the force required to move the leg though the prescribed motion was minimized. The forty-five muscle LifeMOD leg muscle set was then added to both the right and left legs based on relative positions of the hip, knee and ankle. The default attachments of the left quadriceps muscles were modified to insert on the patella based on experimental cadaver measurements (Fig. ).

The measured motion from the dual-limb squat trial was then used in an inverse dynamics simulation to move the lower body as constrained by the hips, ankles, and anatomical knees. During this step the length of each muscle, through their respective via points, was recorded. Next, the motion constraints were removed and each foot was connected to ground with a 6-axis spring-damper. During forward dynamics simulations, muscle forces were generated (within physiological constraints) via feedback control loops that reproduced the muscle lengths recorded during the inverse dynamics training.

Muscle driven forward dynamics squat simulations were run on three versions of the model (Table ). Joint angles, contact forces, muscles forces, meniscus displacement, and ground-reaction forces were recorded during simulation of the three models. Contact pressures on the left knee tibia plateaus were estimated by dividing the contact force on each element by its contact surface area (typically ~ 17 mm2).

RESULTS

Model Validation

Fig. () provides the measured and predicted ground reaction forces during the forward simulation for both the right and left foot for the intact model. The predicted ground reaction forces did not significantly change for the partial and full model versions. Fig. () shows the predicted muscle tension and measured EMG (raw signal rectified and low pass filtered) for four muscles on the left leg where EMG was collected (Peroneus Longus, Gluteus Medius, Vastus Lateralis, Medial Gastrocnemius). Between model versions the magnitude of muscle tension changes, but the overall pattern is generally consistent.

The tibio-femoral compressive force at flexion was 830 N on the lateral side and 2020 N on the medial side at flexion. The total compressive force at the joint was 2850 N or 4.5 times body weight of the subject used in the gait measurements.

Model Comparisons

During the squat cycle, all three model versions maintained the same hip flexion angle (Fig. ). Hip flexion angle varied from 8° at the start of the squat (0.25 s) to 50° at the deepest portion of the squat (1.625 s) and the knee flexion angle varied from approximately 16° to 87° during the same time frame. Table provides the peak predicted contact pressures on the medial and lateral plateaus at extension and flexion for all three model versions. At flexion (87°), peak contact pressure increased on the medial plateau with partial and full ACL transection. On the lateral side, the peak cartilage contact pressure occurred underneath the meniscus at flexion for both the partially and fully transected ACL. At extension, the peak contact force occurred under the lateral meniscus for the full ACL transection.

On the medial plateau, the location of the peak contact pressure and the contact location moved posterior and lateral with partial and full ACL transection at extension (Fig. and Fig. ). At flexion, the location of the peak contact force did not move, but the magnitude of the peak contact force increased. On the lateral plateau, the peak contact location moved to underneath the meniscus at extension for the full transection (Fig. and Fig. ). At flexion, the location of peak contact pressure moved from uncovered cartilage to underneath the meniscus for both the partial and full transection cases (Fig. and Fig. ).

Fig. () provides the motion of the center lateral meniscus rigid body relative to the tibia plateau. The meniscus moves in the posterior direction during extension prior to beginning the squat for all three models. The posterior motion was 3.5 mm greater for both the partial and full transection cases than for the intact knee. The lateral meniscus was pulled in the posterior direction by movement of the femur. The posterior motion occurred faster for the full transection compared to the partial transection case. For all three models, the menisci start in the same position at time zero. There was also a 1 mm difference in medial motion for the transected cases compared to the intact knee. The medial motion was created by tension in the posterior horn attachments and deformation of the entire meniscus as it was stretched in the anterior-posterior direction. Fig. () shows the tibial plateau for all three models at 0.25 second intervals during the squat cycle. At 0.50 seconds the greater posterior motion of the lateral tibia can be seen for the transected cases. In Fig. () the internal force experienced by the menisci rigid bodies is shown at 0.25 second intervals during the motion. After 0.75 seconds the internal forces experienced by the lateral meniscus is altered for the partial and full transection cases.

DISCUSSION

During flexion, the lateral meniscus becomes “wedged” between the tibia and femur for both partial and full ACL transection as indicated by the cartilage contact pressure on elements under the menisci and the changes in the forces experienced by the menisci rigid bodies. This “wedging” effect does not occur for the intact knee. The location of the highest “wedging” force corresponds to the periphery of the central-anterior region of the lateral meniscus. The wedging also occurred at extension for the full transection case at the interior of the central-posterior region of the lateral meniscus (Fig. and Fig. ). The internal forces in the lateral meniscus are also altered for the partial and full transection cases. In particular, the circumferential forces decrease in the posterior lateral meniscus as the knee flexes (Fig. ) and the internal shearing forces increase in the anterior lateral meniscus, where the wedging occurs.

A study evaluating the location of meniscus tears on 476 patients with ACL injury showed that tears were equally distributed between the lateral and medial side and more tears occurred on the posterior of the lateral meniscus than the anterior [20]. The location of meniscus injury evaluated on 159 patients within 3 days of ACL injury during skiing revealed that 83 % of the tears occurred on the lateral side [21]. In the current study, only a simple dual limb squat was simulated and contact pressures were greatest during flexion. More complicated activities such as gait and stair climbing would result in different tibio-femoral loading and motion possibly causing the meniscus “wedging” to occur at different locations.

The tibio-femoral compressive forces in this model were higher than the ones reported by Mundermann et al., [22] for a similar squat motion. In Mundermann’s study the compressive force was measured directly using an instrumented tibial insert after knee reconstruction surgery and it was 2.5 times body weight at flexion. The discrepancy between the natural knee model compressive force and the measured prosthetic compressive force may be explained by the fact that the musculoskeletal model was a combination of a generic skeletal model based on gait subject height and weight and an MRI derived anatomical knee. The stitching of the geometries could have resulted in slight misalignments of the muscle lines of actions that in turn caused higher loads on the cartilage contacts. This issue will be investigated further in future models where the geometries of all bones will be derived from one subject’s MRI data.

Peak contact pressure and contact locations are similar for the partial tear (posterolateral bundle transection) and complete ACL transection cases at flexion, particularly on the lateral side. For the squat simulations, the tibio-femoral contact location shifts to the posterior and lateral direction with ACL transection. This displacement is only a few millimeters and is larger at extension. The change in peak contact and contact patch location on the tibia plateau is consistent with experimental measurements and observations [9, 10, 23]. For both the partial and full transection simulations, small changes in femur motion relative to the tibia push the lateral meniscus in the posterior direction, exposing it to “wedging”. This “wedging” effect increases the loads on the meniscus and it is possible that the increased loading damages the tissue and increases the risk of cartilage degeneration during ACL deficiency [24, 25].

For all simulations only the lower body was modeled. As seen in Fig. (), the upper body of the test subject leans forward during the squat, allowing the body’s center of pressure to stay between its base of support. During the forward dynamics simulations, a LifeMOD tracking agent was used to maintain balance during the squat. The tracking agent is a 6-axis spring located between a “dummy” rigid body and the pelvis. The “dummy” rigid body is driven by a motion constraint that follows pelvis motion measured during the inverse dynamics simulation. If the forward dynamics pelvis motion closely follows the inverse dynamics pelvis motion, the tracking agent will have minimal influence on the model. Also, the tracking agent imparts no force in the vertical direction, regardless of pelvis motion. Measured ground reaction forces are symmetrical in the medial-lateral and anterior-posterior directions, but the simulation results include asymmetry in both axes (Fig. ). For example, both the right and left knee have posterior shear forces while the knee was flexing and an anterior ground reaction shear force when the knee was extending. The model pushes against the tracking agent when these asymmetries occur and these forces will have an influence on muscle force prediction. Misalignment of the knee relative to the hip and ankle and slight misalignment between the left and right sides may be responsible for the predicted ground reaction force asymmetries. Previous work that used a hinge joint for the right knee showed significant asymmetries in the vertical ground reaction forces during the squat [26]. Changing the right knee from a hinge joint to an anatomical joint eliminated the vertical asymmetry. In addition, the vertical ground reaction force predictions may be further improved by including representation of the upper body. In the current model simulations, vertical ground reaction forces were greater for the left foot than the right. Experimentally, the right foot had a greater vertical reaction force (Fig. ). Ongoing work in our lab has shown that including representation of the upper body improves prediction of ground reaction forces, particularly the vertical force distribution between the right and left leg. Future musculoskeletal models will include representation of the upper body.

Combining a cadaver knee with in vivo gait measurements from a subject of similar height and weight allows in vitro validation of the knee model, but it is also a limitation of the study. Additionally, the same measured motion from a subject with an intact ACL was used for the inverse dynamics solution of all three model versions. ACL deficient subjects demonstrate altered gait patterns [27] and it is possible that the kinematics of the dual limb squat would also be different. This altered motion could affect simulation predictions at the knee. Finally, only a limited number of muscle activations were measured with EMG and maximum voluntary muscle contraction was not measured.

The primary objective of this study was to create a musculoskeletal model of the lower extremities that included anatomical representation of the knees, validate the model with experimental data, and then simulate two ACL deficient conditions for a simple motion (squat). The current study is a first step in the development of anatomically correct subject specific musculoskeletal models of the lower extremities to study joint function. Future refinements of this modeling method will include finer segmentations of the menisci and cartilage structures and validation of the models in other dynamic motions.

Supplementary material is available on the publisher’s web site along with the published article.

Table 1.

Stiffness Matrix Parameters for the Lateral and Medial Menisci

Stiffness Parameter	Lateral Meniscus	Medial Meniscus
Kθ (N/mm)	630	540
Kr (N/mm)	592	540
Kz (N/mm)	592	540
Kθr (N/mm)	54	96
Kθz (N/mm)	54	96
Krz (N/mm)	54	96
Tθ (Nmm/deg)	269	222
Tr (Nmm/deg)	378	360
Tz (Nmm/deg)	378	360

Table 2.

Model Versions

Model Name	Procedure
Intact	Intact knee
Partial	Posterolateral bundle of the left knee ACL removed
Full	Anteromedial and posterolateral bundles of the left knee ACL removed

Table 3.

Peak Contact Pressures (MPa) on the Tibia Plateaus

Model Version	Medial		Lateral
	Extension	Flexion	Extension	Flexion
Intact	7.8	14.3	8.5	9.8
Partial	6.4	14.9	8.4	16.4*
Full	7.8	15.3	8.3*	15.5*

Location of peak cartilage pressure under meniscus