Medical Engineering & Physics 36 (2014) 1122–1133

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

Subject-specific evaluation of patellofemoral joint biomechanics during functional activity Massoud Akbarshahi a , Justin W. Fernandez b,c,∗ , Anthony G. Schache a , Marcus G. Pandy a a b c

Department of Mechanical Engineering, University of Melbourne, Victoria 3010, Australia Auckland Bioengineering Institute, University of Auckland, Auckland 1010, New Zealand Department of Engineering Science, University of Auckland, Auckland 1010, New Zealand

a r t i c l e

i n f o

Article history: Received 13 August 2012 Received in revised form 3 June 2014 Accepted 13 June 2014 Keywords: Knee-joint loading Cartilage stress Contact pressure Osteoarthritis X-ray fluoroscopy

a b s t r a c t Patellofemoral joint pain is a common problem experienced by active adults. However, relatively little is known about patellofemoral joint load and its distribution across the medial and lateral facets of the patella. In this study, biomechanical experiments and computational modeling were used to study patellofemoral contact mechanics in four healthy adults during stair ambulation. Subject-specific anatomical and gait data were recorded using magnetic resonance imaging, dynamic X-ray fluoroscopy, video motion capture, and multiple force platforms. From these data, in vivo tibiofemoral joint kinematics and knee muscle forces were computed and then applied to a deformable finite-element model of the patellofemoral joint. The contact force acting on the lateral facet of the patella was 4–6 times higher than that acting on the medial facet. The peak average patellofemoral contact stresses were 8.2 ± 1.0 MPa and 5.9 ± 1.3 MPa for the lateral and medial patellar facets, respectively. Peak normal compressive stress and peak octahedral shear stress occurred near toe-off of the contralateral leg and were higher on the lateral facet than the medial facet; furthermore, the peak compressive stress (11.5 ± 3.0 MPa) was higher than the peak octahedral shear stress (5.2 ± 0.9 MPa). The dominant stress pattern on the lateral patellar facet corresponded well to the location of maximum cartilage thickness. Higher loading of the lateral facet is also consistent with the clinical observation that the lateral compartment of the patellofemoral joint is more prone to osteoarthritis than the medial compartment. Predicted cartilage contact stress maps near contralateral toe-off showed three distinctly different patterns: peak stresses located on the lateral patellar facet; peak stresses located centrally between the medial and lateral patellar facets; and peak stresses located superiorly on both the medial and lateral patellar facets. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Patellofemoral joint (PFJ) pain is a common problem experienced by active adolescents and adults, with incidence rates varying from 9% to 15% [1,2]. The condition is characterized by a gradual onset of peri-patellar pain and these symptoms are often aggravated by daily activities such as walking and stair ambulation [3]. In severe cases, PFJ pain can limit physical activity or result in its cessation altogether. Although the precise etiology of PFJ pain is currently unknown, high joint-contact stress, patellofemoral malalignment, quadriceps muscle weakness, and delayed onset of vastus medialis muscle activity are all thought to be contributing

∗ Corresponding author at: Auckland Bioengineering Institute, UniServices House, University of Auckland, 70 Symonds Street, Auckland 1010, New Zealand. Tel.: +64 9 373 7599. E-mail address: [email protected] (J.W. Fernandez). http://dx.doi.org/10.1016/j.medengphy.2014.06.009 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

factors [4,5]. Importantly, it has been suggested that PFJ pain may be a precursor to the development of PFJ osteoarthritis (OA) [6], which is a common form of OA in the knee joint [7]. Mechanical loading of the patellar cartilage is believed to play a major role in the initiation and progression of PFJ pain symptoms. In vitro experiments [8,9] as well as rigid-body and deformablebody [10] computational models have been used to investigate PFJ biomechanics. In vitro experiments have provided useful insights, but these studies are limited by the fact that they do not replicate the levels of muscle and joint loading observed during daily activity. Similarly, rigid-body models are limited to the study of PFJ kinematics and contact forces. By comparison, deformablebody contact models are capable of calculating cartilage stress distribution, which is likely to be important in identifying factors contributing to structural joint deterioration [11]. Studies aimed at computing PFJ contact stress in vivo are scarce. Bretcher et al. [12] used a model that combined in vivo kinematic and kinetic data with contact areas calculated from magnetic

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

resonance (MR) images and quadriceps moment arms obtained from the literature to estimate PFJ contact forces and stresses during stair ascent. They reported peak PFJ compressive forces and stresses on the order of 3.5 times body weight (BW) and 7 MPa, respectively. These results were derived using a simple model based on regression equations to estimate quadriceps muscle forces but contact stress maps were not computed. Chinkulprasert et al. [13] also used a combination of in vivo joint kinematic and kinetic data, muscle forces obtained from a regression model, and contact areas obtained from the literature to estimate PFJ forces and stresses during forward and lateral step-up and step-down tasks. Peak compressive forces and stresses for a forward step-down task were found to be 51.1 ± 2.7 N/kg and 13.8 ± 0.4 MPa, respectively. Besier et al. [14] developed a subject-specific deformable-body model of the PFJ to estimate patellar cartilage stresses in vivo. An open magnetic resonance (MR) imaging machine was used to measure joint kinematics, whilst muscle forces were estimated using an EMG-driven model. Farrokhi at al. [15] followed a similar approach to calculate joint contact stresses in both healthy subjects and patients with PFJ pain. A major limitation of each of these studies is that the kinematic and kinetic measurements were acquired at different times. The purpose of the present study was to integrate existing capabilities in computational modeling with simultaneous measurements of in vivo joint kinematics and kinetics to estimate subject-specific PFJ kinematics, contact stresses, and stress distribution maps during weight-bearing activity. Stair ascent was investigated because this task is associated with higher PFJ loads than those present during either level walking or stair descent [12,16]. Our specific aims were firstly, to estimate the contact forces and stress distributions on the medial and lateral facets of the patellar cartilage during stair ascent; and secondly, to investigate the correlation between cartilage thickness maps and contact stress maps. We hypothesized that the locations of peak contact stress computed for the PFJ would correspond with the measured regions of maximum patellar cartilage thickness. 2. Methods 2.1. Participants Four healthy adult males (age, 30.5 ± 3 yrs; weight, 71.3 ± 7 kg; height, 178 ± 2 cm) with no history of lower-limb injury gave their informed consent to participate in the study after approval was obtained from the Human Research Ethics Committee at The University of Melbourne. 2.2. Magnetic resonance imaging MR images of each subject’s left lower-limb were acquired using two alternative MR sequences to obtain all the necessary anatomical information. The lower limb from the hip to the ankle was imaged using a T2 fat-suppressed sequence (TE = 12 ms, TR = 23 ms, NEX = 1, slice thickness = 1 mm × 1 mm × 1 mm). These images provided the bony geometries of the femur, tibia, fibula, and patella, as well as the origin and insertion locations of the quadriceps muscles. The patellar cartilage was separated into medial and lateral compartments using the patellar median ridge. A T2 fat-suppressed sequence (TE = 16 ms, TR = 29 ms, NEX = 1, resolution = 0.5 mm × 0.44 mm × 0.44 mm, FOV = 15 cm × 15 cm) using a knee coil was then obtained to visualize the articular cartilage. The MR images were segmented using a thresholding algorithm available in 3D Doctor (Able Software Corp, MA, USA) together with manual intervention. The geometry of the articular cartilage was segmented in two separate layers: the subchondral layer and the

1123

Fig. 1. Top panel: Magnetic resonance (MR) image illustrating how articular cartilage at the tibiofemoral and patellofemoral joints was segmented into two separate layers: a subchondral layer (SL) and an articular layer (AL). Bottom panel: Once the two layers were segmented from the MR images, surface meshes were wrapped to the segmented points on the femoral condyles and patellar facets.

articular layer (Fig. 1). The surface geometries of the bones and articular cartilage were created from the segmented data points. The articular cartilage surfaces were registered to their corresponding subchondral bones using a rigid-body transformation that minimized the Euclidian distances between the surfaces of the bone and the subchondral layer of the articular cartilage. Patellar cartilage thickness maps were generated by calculating the normal distance in the anterior–posterior direction between the two articular cartilage layers using Geomagic (Research Triangle Park, NC). MR-based anatomical coordinate systems for the femur, tibia, and patella were defined as described in Fernandez et al. [17]. 2.3. Gait experiments Each subject performed a stair-ascent task on a customdesigned staircase comprised of four steps (step height = 31 cm; step width = 79 cm; step depth = 40 cm). Seventeen reflective markers were mounted on the subject’s pelvis, left thigh, left shank, and left foot. A video motion capture system (VICON, Oxford Metrics Inc.) with nine cameras sampling at 120 Hz was used to record the marker-based kinematic data. Ground reaction forces were measured using two portable force platforms (AMTI Accugait, AMTI Corporation, Watertown, MA) mounted on the stairs. Muscle EMG data were recorded using a telemetry system (BIOTEL99, Neurodata, Vienna, Austria). Surface electrodes (Kendall Medi-trace 100, Tyco Healthcare Group) were placed over the rectus femoris, vastus medialis, and vastus lateralis muscles of the subject’s left leg. EMG onset and offset times were determined by applying a

1124

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

Fig. 2. Diagram illustrating how gait analysis experiments, X-ray fluoroscopy, and MR-based subject-specific modeling were used to calculate lower-limb muscle forces and estimate patellofemoral joint loading in vivo.

Teager–Kaiser energy (TKE) filter to the raw EMG signal [18,19]. This method has previously been used to quantify muscle EMG onset and offset times during walking and running [20]. Two synchronized X-ray units were positioned in series to capture dynamic X-ray images of the knee joint at 30 Hz (Fig. 2). Each X-ray system was calibrated and the X-ray images corrected for distortion. Tibiofemoral joint (TFJ) kinematics were measured by registering the 3D bone models to the 2D X-ray images using a feature-based pose-estimation algorithm described previously [21]. The accuracy of the motion capture data has been previously reported to range in magnitude from 4◦ to 13◦ in rotation and 13.9 to 16.1 mm in translations [22]. The accuracy of the fluoroscopic measurements was estimated to be 0.6 degrees for rotations, 1.5 mm for anterior–posterior translations and 3 mm for medial–lateral translations [21]. 2.4. Musculoskeletal model of the body A three-dimensional muscle-actuated model of the body was used to calculate lower-limb muscle forces for one cycle of stair ascent [23]. The skeleton was represented as an 8-segment, 21-degree-of-freedom linkage. The head, arms and torso were modeled as a single rigid body, which articulated with the pelvis via a ball-and-socket back joint. Each hip was modeled as a ball-andsocket joint, each knee as a hinge, each ankle-subtalar complex as a universal joint, and each metatarsal joint as a hinge. The wholebody model was actuated by 58 Hill-type muscle-tendon units. Subject-specific musculoskeletal models were created by scaling the segmental inertial properties of the model, muscle attachment sites, and muscle paths to each subject’s height and weight. Specifically, virtual markers were defined in the generic musculoskeletal model based on the positions of the reflective markers attached to each subject. The anthropometric properties of the body segments and muscle-tendon units were then linearly scaled using the relative distances between the virtual markers and reflective markers. The following body-segment and muscle-tendon properties were

scaled: the length, width, and depth of each body segment; location of the center of mass and mass moment of inertia of each body segment; origin and insertion sites of each muscle-tendon unit; and the locations of muscle-tendon wrapping over bone. 2.5. Muscle force calculations Muscle forces were found using inverse dynamics and static optimization. Measurements of the subject’s motion and ground reaction forces were input into the corresponding subject-specific model, and inverse dynamics was used to calculate the net moments exerted about the back, hip, knee and ankle joints for one stair-ascent cycle. The net joint moments were decomposed into individual muscle forces by solving an optimization problem that minimized the sum of the squares of muscle activations. The optimization solution was constrained to the force–length–velocity surface of each muscle [24]. As a means of verifying the model calculations, the sequence and timing of muscle forces predicted by the model were compared against EMG measurements obtained for each subject. This approach has been used in a number of previous studies that have calculated lower-limb muscle forces during human locomotion; see Pandy and Andriacchi [25] for review. For each muscle the percentage of predicted force that correlated temporally with EMG onset time was calculated to quantify the agreement between measured muscle EMG activity and modelpredicted muscle force. Similarly, the percentage of EMG onset time that did not correlate temporally with the calculated muscle force was used to quantify the lack of agreement between measured muscle EMG activity and model-predicted muscle force. 2.6. Finite-element model of the patellofemoral joint A deformable finite-element (FE) model was developed to calculate PFJ contact forces and stresses. Only contact between the femur and patella was modeled. Uniform hexahedral element meshes of the femur, tibia and patella bones as well as the articular

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

1125

Fig. 3. Diagram showing the muscles actuating the finite-element model of the patellofemoral joint developed in this study. Four quadriceps muscles actuated the model: vastus medialis, VM; vastus lateralis, VL; vastus intermedius, VI; and rectus femoris, RF. The patellar tendon was represented by three bundles. Four bundles were used to model the actions of the medial and lateral retinaculum: medial retinaculum bundle, MRB; lateral retinaculum bundle, LRB. The inset shows how the different insertion locations of the patellar tendon were altered in the model. M, medial shift; L, lateral shift; P, proximal shift; and D, distal shift.

cartilages of the distal femur and patella were created using TrueGrid (XYZ Scientific Ltd., Livermore, CA). Articular cartilage was modeled as a deformable, linear, elastic material (modulus of elasticity, E = 10 MPa; Poisson’s ratio,  = 0.45) [9,26]. The bones were modeled as rigid shells for computational efficiency, as the focus of the present study was to estimate cartilage deformation. The resultant vasti force obtained from the whole-body model was partitioned across the three sub-regions of this muscle group – vastus lateralis, vastus intermedius, and vastus medialis – based on measurements of the cross-sectional areas of these sub-regions obtained from MRI. The cross-sectional area of each muscle was computed by dividing the volume of the muscle by its arc-length. Muscle lines of action were defined by computing the centroidal paths of the muscles, which was reproduced directly from the MR images acquired for each subject and represented in the FE model as connector elements (Fig. 3). The superior end of each muscle was rigidly fixed to the femur at its attachment point, which was determined from the MR images. Contact between each muscle’s line-of-action and the femur was modeled to account for wrapping of the muscle-tendon unit around the bone in deep flexion. The patellar tendon was represented as three separate nonlinear bands with the force–elongation curve taken from cadaveric data [27]. To estimate the length of the patellar tendon, the sagittal-plane position of the patella was matched to that measured from the fluoroscopic images at heel-strike (i.e., the beginning of the simulated stair-ascent cycle). Four additional nonlinear bands were attached to the medial and lateral sides of the patella and femoral condyles to model the effects of skin, retinaculum, and the surrounding soft tissues in the anterior compartment of the knee, including the knee capsule. The mechanical properties of the knee capsule and retinaculum were adapted from the model reported by Shelburne and Pandy [28].

The FE model was developed in Abaqus (Simulia, Dassault Systèmes, Providence, RI) as a quasi-static explicit simulation with only force and moment equilibrium considered at each time point. The model simulated contact between the patella and the femur only. An augmented Lagrangian contact method was chosen with a ‘hard-contact’ constraint to ensure there was no cartilage overlap for PFJ contact. The simulation was executed at 60 discrete time points with a fixed femur and moving tibia. The inputs to the model were the active muscle forces obtained from the whole-body musculoskeletal model, the passive forces applied by the patellar tendon and the soft tissues attached to the patella, and the six degrees of freedom of tibiofemoral joint motion measured from X-ray fluoroscopy. The model outputs included the medial and lateral patellar contact forces, normal (compressive) and octahedral shear stresses, and PFJ kinematics (i.e., all six degrees of freedom defining patellar position). PFJ kinematics was predicted from the equilibrium of active muscle forces, passive soft-tissue forces, and contact forces. To partially validate the PFJ kinematics predicted by the model, the PFJ flexion angle was also measured from fluoroscopy and compared to the corresponding angle calculated in the model. All kinematics were calculated using the Joint Coordinate System proposed by Grood and Suntay [29]. We evaluated the sensitivity of the FE model calculations to changes in the input muscle forces and model parameters, specifically, quadriceps forces, patellar-tendon insertion locations, and the elastic modulus of articular cartilage. The forces applied by the vastus medialis, vastus lateralis, vastus intermedius and rectus femoris were each increased by 10% in turn, while the remaining muscle forces were adjusted to keep the magnitude of the net knee-joint moment equal to that calculated from inverse dynamics. Patellar tendon insertion locations were moved ±5.5 mm inferosuperiorly and ±2.6 mm mediolaterally to simulate the error involved

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

HS

10 5 0 -5 -10 -15

0

20 40 60 80 Stair Rise Cycle % Gait Cycle (%)

100

HS CTO 15

TO

HS

Medial Tibial Translation (mm)

TO

Superior Tibial Translation (mm)

HS CTO 15

10 5 0 -5 -10 -15

HS CTO 15

0 -5

-10

0

20 60 40 80 StaiGait r Rise Cycle Cycle (%)%

-15

100

100

5

5

40 20 0 -20 0

Internal Rotation (Deg)

10

60

0 -5 -10 -15 -20

100

0

20 40 60 80 StaGait ir Rise Cycle Cycle (%)%

100

0

20 40 60 80 StairGait RisCycl e Cycle e (%)%

100

0 -5 -10 -15 -20 -25

-25 20 40 60 80 StairGait Rise Cycle % Cycl e (%)

HS

5

10

80

TO

10

120

Abduction (Deg)

Flexion (Deg)

Anterior Tibial Translation (mm)

1126

0

20 40 60 80 StairGait Rise CyCycle cle (%)%

100

SUBJECT1 SUBJECT2 SUBJECT3 SUBJECT4 Fig. 4. Tibiofemoral translations and rotations measured from X-ray fluoroscopy plotted as a percentage of the stair-ascent cycle for each of the four subjects. Positive displacements represent anterior, proximal, and medial translations of the tibia relative to the femur, while negative displacements represent posterior, distal, and lateral translations of the tibia relative to the femur. Positive rotations represent flexion, abduction, and internal rotation of the tibia relative to the femur, while negative rotations represent extension, adduction, and external rotation of the tibia relative to the femur. Prominent events of the stair-ascent cycle indicated on the graphs are as follows: HS, heel-strike of the ipsilateral leg; CTO, contralateral toe-off; and TO, toe-off of the ipsilateral leg.

in identifying tendon attachment sites from MRI. The elastic modulus of cartilage was perturbed by ±10% and ±20% from its nominal value (E = 10 MPa). FE model convergence was also investigated using four different mesh densities: 3.0 mm, 2.25 mm, 1.7 mm and 1.5 mm. The model sensitivity calculations were performed on one subject-specific model and the results for peak contact force, peak contact stress, and contact area compared at contralateral toe-off. 3. Results A 2.25 mm mesh provided sufficient resolution for determining all stress measures as well as the resultant contact force acting at the PFJ (Table 1). Contact forces were minimally affected by a change in the modulus of elasticity assumed for articular cartilage. However, contact areas and stresses changed by a maximum of 10% and 7% when the elastic modulus was altered by 10% and 20%, respectively (Table 1). Contact areas and forces were sensitive to the shift in the location of the patellar tendon on the tibial tuberosity; for example, a 2.6 mm shift in the mediolateral direction changed the peak medial contact force by 19% and the medial contact area by 15%. Patellofemoral contact force and contact stress were less sensitive to changes in rectus femoris and vastus intermedius forces than to variations in the forces applied by vastus lateralis and vastus medialis. Overall, the peak PF contact force and compressive contact stress at contralateral toe-off increased by less than 0.5%.

However, when vastus lateralis force was increased by 10%, the force acting on the lateral patellar facet increased by 2.9% while that acting on the medial patellar facet decreased by 3.2%; the corresponding total peak compressive stress increased by 1.1%. Similarly, when vastus medialis force was increased by 10%, the force acting on the medial patellar facet increased by 4% whereas that acting on the lateral patellar facet decreased by 1.7%; the total peak compressive stress decreased by 1.1%. Subjects exhibited similar patterns of sagittal-plane TFJ motion (derived from fluoroscopic data) albeit with differences in the magnitudes of the peak joint translations (Fig. 4). In general, the femur translated anteriorly and superiorly with knee extension, and posteriorly and inferiorly with knee flexion. There was no consistent pattern of motion for either TFJ abduction–adduction or internal–external rotation. The range of motion for abduction–adduction was less than 7◦ for all subjects. Trends in the model-predicted PFJ kinematics were also reasonably consistent across subjects (Fig. 5). The model patella flexed and translated anterosuperiorly with knee extension. The non-sagittalplane rotations of the patella (tilt and twist) displayed amplitudes of less than 15◦ . No consistent patterns of motion were observed for patellar tilt and patellar twist as a function of knee flexion. The average root-mean-square (RMS) difference between patellar flexion angle predicted by the model and that measured from the fluoroscopy images was 2.8◦ across all subjects. The average RMS errors ranged from 0.4◦ (subject 1) to 4.1◦ (subject 3) (Fig. 6).

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

HS

40

20

0 0

20

40 80 60 Stair Rise Cycle % Gait Cycle (%)

20

0

-20

HS CTO 40

-40

100

0

40

HS

TO

HS CTO 40

Medial Patellar translation (mm)

TO

Proximal Patellar Translation (mm)

Anterior Patellar Translation (mm)

HS CTO 60

20

60 80 40 Stair Rise Cycle % Gait Cycle (%)

1127

HS

20

0

-20

-40 0

100

TO

20

20

10

10

20

40 60 80 100 Stair Rise Cycle % Gait Cycle (%)

0 -20 -40

Patellar Twist (deg)

Patellar Tilt (deg)

Patellar Flexion (deg)

20

0

-10

0

-10

-60 -80 0

-20 20

40

60

80

100

0

20

40

60

80

100

-20

0

Stair GaiRise t CyclCycle e (%) %

Stair GaiRise t CycCycle le (%) %

20

40

60

80

100

Stair GaitRi Cse e (%) % yclCycle

SUBJECT1 SUBJECT2 SUBJECT3 SUBJECT4 Fig. 5. Patellofemoral translations and rotations predicted by the finite-element model plotted as a percentage of the stair-ascent cycle for each of the four subjects. Positive displacements represent anterior, proximal, and medial translations of the patella relative to the femur, while negative displacements represent posterior, distal, and lateral translations of the patella. Positive rotations represent flexion, medial tilt, and medial twist of the patella relative to the femur, while negative rotations represent extension, lateral tilt, and lateral twist of the patella relative to the femur. Prominent events of the stair-ascent cycle indicated on the graphs are as follows: HS, heel-strike of the ipsilateral leg; CTO, contralateral toe-off; and TO, toe-off of the ipsilateral leg.

The force calculated for vasti was similar in shape across all subjects, and the timing of predicted muscle activity was consistent with EMG measurements obtained from experiment. Specifically, the peak force developed by the vasti occurred near contralateral

toe-off for all subjects (Fig. 7), with the force developed by vastus lateralis highest and that developed by rectus femoris force lowest (2.7 ± 0.5 BW and 0.9 ± 0.1 BW, respectively). Across all four subjects the percentage of predicted muscle force that temporally

Table 1 Effect of mesh size on FE model convergence and the sensitivity of changes in the values assumed for patellar tendon insertion locations, elastic modulus of articular cartilage, and quadriceps muscle force on FE model calculations. Data shown in bold typeface represent the nominal condition assumed in the model. Convergence/sensitivity analysis

Peak stress (MPa)

Peak force (N)

Area (mm2 )

CTO

Lat

Med

Total

Lat

Med

Total

Change in average mesh size (mm)

3.0 × 3.0 Nominal 2.25 × 2.25 1.7 × 1.7 1.5 × 1.5

13.2 14.4 14.6 14.8

5477 5059 4968 4939

411 560 598 648

5926 5661 5566 5587

511 465 457 451

176 160 158 156

687 625 615 607

Change in modulus of elasticity

Nominal −20% −10% +10% +20%

14.4 13.4 13.9 14.9 15.3

5059 5024 5043 5077 5095

560 586 571 566 548

5661 5610 5613 5644 5643

511 540 521 497 472

176 179 176 160 158

687 720 696 657 630

Change in patellar tendon insertion location

Nominal Medial Lateral Proximal Distal

14.4 14.0 14.8 14.7 14.1

5059 4872 5250 4944 5172

560 669 511 653 561

5661 5541 5761 5597 5733

511 499 511 499 509

176 183 148 176 170

687 683 659 677 679

10% increase in quadriceps muscle force

Nominal VL VM VI RF

14.4 14.56 14.24 14.46 14.43

5059 5206 4973 5079 5069

560 542 582 562 561

5661 5830 5654 5684 5672

511 526 502 513 512

176 170 183 177 176

687 696 685 690 688

1128

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

HS

CTO

TO

HS

Patellar (deg) Patellar Flexion Flexion (deg)

Patellar Flexion (deg)

CTO

TO

30

30 10 -10 -30 -50

10 -10 -30 -50

Subject 2

Subject 1 -70

0

10

20

30

40

-70

50

10

0

Stair Rise Cycle Gait Cycle (%) %

30

40

Patellar Flexion Flexion (deg) Patellar (deg)

10 -10 -30 -50

10 -10 -30 -50

Subject 3 -70

50

30

30

Patellar Flexion (deg)

20

Stair Rise Cycle Gait Cycle (%) %

0

10

20

30

40

Subject 4 -70

50

0

10

Gait Rise CycleCycle (% ) % Stair

20

30

40

50

GaRise it CycCycle le (%) % Stair

Fluoroscopy Kinematics FE Model Fig. 6. Patellofemoral flexion angle measured from dynamic fluoroscopy and compared to that predicted by subject-specific FE models of the PFJ. Results are shown for each of the four participants in this study. Prominent events of the stair-ascent cycle indicated on the graphs are as follows: HS, heel-strike of the ipsilateral leg; CTO, contralateral toe-off; and TO, toe-off of the ipsilateral leg.

Table 2 Comparison of predicted muscle force and measured EMG activity for all four subjects. The first four columns show the percentage of predicted muscle force that overlaps temporally with measured EMG activity for the rectus femoris (RF), vastus medialis (VM) and vastus lateralis (VL) muscles for each subject. The second four columns show the percentage of measured EMG activity that does not overlap with predicted muscle force for the same three muscles. % of predicted muscle force that overlaps with EMG activity

RF VM VL

% of EMG activity that does not overlap with any predicted muscle force

Subject 1

Subject 2

Subject 3

Subject 4

Subject 1

Subject 2

Subject 3

Subject 4

89 100 74

44 78 91

48 64 100

33 81 95

6 23 0

0 26 5

14 11 27

14 0 5

overlapped with measured EMG activity ranged from 33% to 89%, 64% to 100% and 74% to 100% for the rectus femoris, vastus medialis and vastus lateralis, respectively (Table 2). Further, the percentage of measured EMG activity that did not temporally overlap with predicted muscle force ranged from 6% to 14%, 0% to 26% and 0% to 27%, respectively, for the rectus femoris, vastus medialis and vastus lateralis. Overall, the muscle forces calculated for Subject 1 displayed the closest agreement with EMG measurements, whereas the timing predicted for vastus lateralis (the largest quadriceps muscle) showed the best agreement with experiment. Contact forces were higher on the lateral facet of the patella for all subjects. Peak PFJ contact forces occurred near contralateral toe-off and ranged from 3.7 BW to 6.1 BW on the lateral facet and from 0.7 BW to 1.3 BW on the medial facet (Fig. 8). Peak PFJ contact stresses averaged across subjects were 8.2 ± 1.0 MPa and 5.9 ± 1.3 MPa for the lateral and medial facets, respectively (Fig. 8). Maximum localized normal compressive stress at the maximum loading point near contralateral toe-off ranged from 8.3 MPa to 14.4 MPa across all subjects (Fig. 9). At this time point, the

subject-specific contact stress maps showed three distinctly different patterns: peak loading located on the lateral patellar facet; peak loading located centrally between the medial and lateral patellar facets; and peak loading located superiorly on both the medial and lateral patellar facets. For subjects 2, 3 and 4, the region of high contact compressive stress predicted by the model was similar in location to the region of maximum cartilage thickness measured from the MR images (Fig. 9). Peak octahedral shear stress was lower than the normal compressive stress; the average peak value of shear stress was 5.2 ± 0.8 MPa compared to an average peak value of 11.4 ± 3.1 MPa for normal stress (Fig. 10). The shear stress spanned a larger area than the normal stress, with the maximum shear stress located centrally and laterally. 4. Discussion A subject-specific modeling framework was used to investigate the biomechanical behavior of the PFJ during stair ambulation.

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

1.5

HS

CTO

TO

3.0

HS

1129

CTO

TO

VM

1.0

Force(BW) (BW) Force

Force (BW)

RF

0.5

1.0

0

EMG

0

2.0

3.0

3.0

2.0

1.0

0

10

20 30 40 GaRise it CycCycl le(%)e % Stair

50

2.0

1.0

0

0

10

20

30

40

50

GaRise it CycCycle le(%) % Stair

EMG

0

VI Force (BW) Force (BW)

Force (BW)

VL

SUBJECT1 SUBJECT2 SUBJECT3 SUBJECT4 Fig. 7. Quadriceps muscle forces calculated for each of the four subjects. Horizontal bars shown below each plot indicate the periods of electromyographic (EMG) activity recorded for each muscle as determined by Teager–Kaiser energy filtering of the raw EMG signal. No EMG data were recorded for vastus intermedius (VI) as this is a deep-lying muscle of the thigh. Other muscle symbols: RF, rectus femoris; VM, vastus medialis; VL, vastus lateralis. Prominent events of the stair-ascent cycle indicated on the graphs are as follows: HS, heel-strike of the ipsilateral leg; CTO, contralateral toe-off; and TO, toe-off of the ipsilateral leg. Only data for the stance phase of the stair-ascent cycle are shown.

Geometric bone and cartilage models were developed for four subjects using MR imaging. In vivo kinematic and kinetic data were recorded simultaneously and applied to a deformable, subjectspecific, FE model of the PFJ to compute patellar kinematics, contact forces, and contact stresses. Contact forces were 4–6 times higher on the lateral facet than on the medial facet, and average contact stresses were also generally higher on the lateral facet (Fig. 8). Maximum localized shear and normal stresses were located on the lateral facet of the patella (Figs. 9 and 10). The high stress patterns evident for the lateral facet are consistent with the higher prevalence of OA in this region of the patella [30], supporting the view that mechanical loading plays an important role in the development and progression of PFJ OA [11]. There are several limitations of this study that must be acknowledged. First, and perhaps most significantly, the time and computational cost associated with generating subject-specific FE models of PFJ contact mechanics limited the analysis to four subjects. Second, the use of single-plane fluoroscopy limited the accuracy of the kinematic measurements of knee-joint motion during stair ascent. The errors involved in measuring out-of-plane bone movements using single-plane fluoroscopy can be as high as 3 mm [21], which compromised our measurements of the mediolateral translations of the tibia relative to the femur, and hence the

predicted out-of-plane movements of the patella, particularly near full extension where patellar motion is not constrained by the shape of the trochlear groove. Third, the knee was modeled as a hinge in the whole-body musculoskeletal model, thus ignoring the effects of adduction–abduction and internal–external rotation torques. As a consequence, quadriceps forces in the model may have been underestimated as these quantities were estimated without considering the balance of forces about all three axes of knee-joint motion. Fourth, the whole-body musculoskeletal model did not explicitly include a patella. However, the function of the patella as a spacer and a lever was incorporated in the model by using measured values of the knee-extensor moment arm obtained from cadaver experiments [31]. In addition, the model calculations of patellar contact forces and pressures accounted for wrapping of the quadriceps muscles around the distal femur which occurs at high flexion angles of the knee. Finally, the estimates obtained for the lower-limb muscle forces and knee-joint contact forces cannot be quantitatively validated because there is currently no non-invasive method available for directly measuring these quantities in healthy subjects. Nonetheless, the muscle force predictions were in good temporal agreement with the patterns of measured EMG activity (Fig. 7). In addition, the musculoskeletal and finite-element modeling methods used in the

1130

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

CTO

6

4

2

0

0

10

20

30

40

HS

CTO

1

0.5

0

10

PF Medial Medial Contact PF ContactStress Stress(MPa) (MPa)

4 2

10

20

30

40

20

30

40

50

Stair Gait RiseCycle(%) Cycle %

6

0

TO

1.5

0

50

8

0

2

StairGait Rise Cycle % Cycle(%)

B 10 PF Lateral Contact Stress (MPa)

TO PF Medial Medial Contact PF ContactForce Force(BW) (BW)

PF Lateral Contact Force (BW)

A 8 HS

50

10 8 6 4 2 0 0

StairGRi le)% aise t CyCyc cle(%

10

20

30

40

50

StairGRi aise t CyCyc cle(le %)% SUBJECT1 SUBJECT2 SUBJECT3 SUBJECT4

Fig. 8. Patellofemoral contact forces (A) and average patellofemoral normal compressive stresses (B) predicted by the finite-element model for each of the four subjects.

present study have been assessed previously by Kim et al. [32]. In that study, model predictions of knee muscle forces for walking at three different speeds (slow, normal and fast) were quantitatively evaluated against in vivo measurements of joint contact loading acquired from an instrumented implant. Calculated and measured tibial contact forces were in good agreement for all three walking speeds. Average root mean square (RMS) errors for the medial, lateral, and total contact forces were, respectively, 0.21 ± 0.06 BW, 0.17 ± 0.05 BW, and 0.27 ± 0.07 BW, calculated over one gait cycle and across all walking trials [32]. However, the influence of the passive collateral ligaments, which may affect the behavior of those muscles with large mediolateral moment arms, was not quantified by Kim et al. [32]. The quadriceps muscles do not have large mediolateral moment arms as they insert onto the patella via the quadriceps tendon; thus, we would not expect the collateral knee ligaments to significantly affect the calculated values of quadriceps force, and hence PFJ contact forces and stresses, obtained in the present study. Numerous in vitro [33–35] and in vivo [12,14] studies have investigated PFJ contact mechanics. These studies report values ranging from 229 to 4000 N for patellar contact force, 1.3 to 13.8 MPa for patellar contact (compressive) stress and 0.6 to 2.5 MPa for octahedral shear stress [12–14,34]. The majority of these estimates are substantially lower than those obtained in the present study; for example, the peak contact force and peak average contact stress acting on the lateral patellar facet were found to be 4922 N (6.1 BW) and 9.2 MPa, respectively. This disparity is most likely attributable

to differences in the experimental design adopted in the various studies. In particular, a relatively large step height was used in the current study, which resulted in larger peak knee flexion angles, and hence larger knee extensor torques and quadriceps forces. A recent in vitro study by Goudakos et al. [34] simulated a stair ascent task using knee flexion angles of 12◦ and 30◦ . Muscle forces of a similar magnitude to those shown in Fig. 7 were applied to the cadaver specimens, resulting in magnitudes of joint contact forces and stresses consistent with our results. The spatial distributions of PFJ normal contact stress and octahedral shear stress calculated in the present study were also comparable to those reported previously [14,15], with regions of high stress concentrated on the central and lateral portions of the patellar cartilage. The maximum RMS difference between the measured and calculated values of patellar flexion angle was 4.1◦ . The predicted PFJ kinematic data are in close agreement with in vitro data reported by Ahmed et al. [8]. These investigators reported approximately 50 mm of translation in both the proximodistal and anteroposterior directions as the knee flexes from full extension to 90◦ , which is similar to the results shown in Fig. 5. The amplitudes of patellar flexion, tilt, and twist calculated in the model are also consistent with experimental data reported by Ahmed et al. [8]. Patellar cartilage thickness maps measured in this study ranged from 0 to 6 mm with an average of 2.8 mm, which is similar to data reported by Cohen et al. [36] (range: 0–6 mm; average: 3.05 mm). The patellar cartilage was generally thicker on the lateral facet. Thus, for three of the four subjects, the location of maximum

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

1131

Fig. 9. Left panels: Localized normal contact stress maps for the patella predicted by the finite-element model for each of the four subjects at the maximal loading point near contralateral toe-off. Right panels: Patellar cartilage thickness maps measured for each of the four subjects.

localized normal stress corresponded well to the location of maximum cartilage thickness (Fig. 9). This result is consistent with the findings of previous studies which show that cartilage responds favorably to joint loading in healthy adults [33]. In contrast, the maximum localized normal stress calculated for subject 1 occurred proximal to the region of maximum cartilage thickness and was caused by a more inferior position of the patella (Fig. 9). One mechanism reported for the onset of OA is loading in regions of thinner cartilage, followed by migration to regions of high load in the later phases of the disease when mechanical load plays a destructive role and accelerates the process of cartilage degeneration [11]. In particular, octahedral shear stress can accelerate the process of endochondral ossification (i.e., the transformation of cartilage into bone) and result in further thinning of the cartilage tissue [37]. The normal and shear stresses predicted for subject 1 (Figs. 9 and 10) highlight the importance of subjectspecific loading patterns, as peak loading was observed to occur

in more superior regions of the patellar facets where cartilage is thinner. The variation in the patterns of patellar contact stress obtained for the four participants in this study (Fig. 9) highlights the role that computational modeling may play in subject-specific diagnosis and treatment of musculoskeletal conditions. Subject 1 showed a pattern of split loading between the lateral and medial patellar facets, with slightly higher peak stresses on the medial facet. Subject 2 showed a centrally located stress pattern near the patellar ridge, although the contact stresses were more focused laterally. Subjects 3 and 4 showed peak stresses localized on the lateral facet of the patella. The distribution of contact stress is influenced by the articular geometry that is unique to each individual and cannot be revealed by considering average stress alone. Thus, whilst PFJ lateral compartment loading is most commonly reported, a subject’s geometry, kinematics and muscle loading all give rise to differences in the distribution of contact stress. Differences in patellar loading

1132

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133

Fig. 10. Left panels: Localized octahedral shear stress maps for the patella predicted by the finite-element model for each of the four subjects at the maximal loading point near contralateral toe-off. Right panels: Patellar cartilage thickness maps measured for each of the four subjects.

patterns are evident from both modeling and experimental studies reported in the literature. The FE model calculations performed by Farrokhi et al. [15] showed that subjects with patellofemoral pain have increased loading on the lateral facet of the patella, whereas biomechanical experiments performed on human cadavers showed increased levels of cartilage wear on the medial patellar facet [38]. To our knowledge, this is the first study to investigate the subject-specific in vivo biomechanical behavior of the PFJ during stair ambulation. Contact force, peak compressive stress, and peak octahedral shear stress all occurred near contralateral toe-off and were higher on the lateral facet than the medial facet. The dominant stress pattern on the lateral facet corresponded well to the location of maximum cartilage thickness. The correlation between patellar cartilage thickness maps and patellar contact stress patterns requires further investigation to gain a more complete understanding of the response of articular cartilage to repetitive joint loading.

Ethical approval Approval was obtained from the Human Research Ethics Committee at The University of Melbourne. Acknowledgements This work was supported in part by Australian Research Council Discovery Grants DP0772838, DP1095366 and DP120101973 and a VESKI Innovation Fellowship to M.G.P. Conflict of interest This study has no competing interests. References [1] Milgrom C, Kerem E, Finestone A, Eldad A, Shlamkovitch N. Patellofemoral pain caused by overactivity. A prospective study of risk factors in infantry recruits. J Bone Joint Surg Am 1991;73-A:1041–3.

M. Akbarshahi et al. / Medical Engineering & Physics 36 (2014) 1122–1133 [2] Schwellnus MP, Jordaan G, Noakes TD. Prevention of common overuse injuries by the use of shock absorbing insoles. A prospective study. Am J Sports Med 1990;18:636–41. [3] Crossley K, Bennell K, Green S, McConnell J. A systematic review of physical interventions for patellofemoral pain syndrome. Clin J Sports Med 2001;11:103–10. [4] Fulkerson JP, Shea KP, Ewing JW. Mechanical basis for patellofemoral pain and cartilage breakdown. In: Ewing JW, editor. Articular cartilage and knee joint function: basic science and arthroscopy. New York: Raven Press; 1990. p. 93–101. [5] Witvrouw E, Lysens R, Bellemans J, Cambier D, Vanderstraeten G. Intrinsic risk factors for the development of anterior knee pain in an athletic population. A two-year prospective study. Am J Sports Med 2000;28:480–9. [6] Wilson DR, McWalter EJ, Johnston JD. The measurement of joint mechanics and their role in osteoarthritis genesis and progression. Rheum Dis Clin North Am 2008;34:605–22. [7] Duncan RC, Hay EM, Saklatvala J, Croft PR. Prevalence of radiographic osteoarthritis—it all depends on your point of view. Rheumatology (Oxford) 2006;45:757–60. [8] Ahmed AM, Duncan NA, Tanzer M. In vitro measurement of the tracking pattern of the human patella. J Biomech Eng 1999;121:222–8. [9] Elias JJ, Wilson DR, Adamson R, Cosgarea AJ. Evaluation of a computational model used to predict the patellofemoral contact pressure distribution. J Biomech 2004;37:295–302. [10] Cohen ZA, Roglic H, Grelsamer RP, Henry JH, Levine WN, Mow VC, et al. Patellofemoral stresses during open and closed kinetic chain exercises. An analysis using computer simulation. Am J Sports Med 2001;29:480–7. [11] Carter DR, Beaupre GS, Wong M, Smith RL, Andriacchi TP, Schurman DJ. The mechanobiology of articular cartilage development and degeneration. Clin Orthop Relat Res 2004;427 Suppl.:S69–77. [12] Brechter JH, Powers CM. Patellofemoral joint stress during stair ascent and descent in persons with and without patellofemoral pain. Gait Posture 2002;16:115–23. [13] Chinkulprasert C, Vachalathiti R, Powers CM. Patellofemoral joint forces and stress during forward step-up, lateral step-up, and forward step-down exercises. J Orthop Sports Phys Ther 2011;41:241–8. [14] Besier TF, Gold GE, Beaupre GS, Delp SL. A modeling framework to estimate patellofemoral joint cartilage stress in vivo. Med Sci Sports Exerc 2005;37:1924–30. [15] Farrokhi S, Keyak JH, Powers CM. Individuals with patellofemoral pain exhibit greater patellofemoral joint stress: a finite element analysis study. Osteoarthritis Cartilage 2011;19:287–94. [16] Matthews LS, Sonstegard DA, Henke JA. Load bearing characteristics of the patello-femoral joint. Acta Orthop 1977;48:511–6. [17] Fernandez JW, Akbarshahi M, Kim HJ, Pandy MG. Integrating modelling, motion capture and X-ray fluoroscopy to investigate patellofemoral function during dynamic activity. Comput Methods Biomech Biomed Eng 2008;11:41–53. [18] Solnik S, Rider P, Steinweg K, DeVita P, Hortobagyi T. Teager–Kaiser energy operator signal conditioning improves EMG onset detection. Eur J Appl Physiol 2010;110:489–98. [19] Li X, Zhou P, Aruin AS. Teager–Kaiser energy operation of surface EMG improves muscle activity onset detection. Ann Biomed Eng 2007;35:1532–8.

1133

[20] Dorn TW, Schache AG, Pandy MG. Muscular strategy shift in human running: dependence of running speed on hip and ankle muscle performance. J Exp Biol 2012;215:1944–56. [21] Akbarshahi M, Schache AG, Fernandez JW, Baker R, Banks S, Pandy MG. Noninvasive assessment of soft-tissue artifact and its effect on knee joint kinematics during functional activity. J Biomech 2010;43:1292–301. [22] Benoit DL, Ramsey DK, Lamontagne M, Xu L, Wretenberg P, Renstrom P. Effect of skin movement artifact on knee kinematics during gait and cutting motions measured in vivo. Gait Posture 2006;24:152–64. [23] Anderson FC, Pandy MG. A dynamic optimization solution for vertical jumping in three dimensions. Comput Methods Biomech Biomed Eng 1999;2: 201–31. [24] Anderson FC, Pandy MG. Static and dynamic optimization solutions for gait are practically equivalent. J Biomech 2001;34:153–61. [25] Pandy MG, Andriacchi TP. Muscle and joint function in human locomotion. Annu Rev Biomed Eng 2010;12:401–33. [26] Kim YJ, Bonassar LJ, Grodzinsky AJ. The role of cartilage streaming potential, fluid flow and pressure in the stimulation of chondrocyte biosynthesis during dynamic compression. J Biomech 1995;28:1055–66. [27] Reeves ND, Maganaris CN, Narici MV. Effect of strength training on human patella tendon mechanical properties of older individuals. J Physiol 2003;548:971–81. [28] Shelburne KB, Pandy MG. A musculoskeletal model of the knee for evaluating ligament forces during isometric contractions. J Biomech 1997;30:163–76. [29] Grood ES, Suntay WJ. A joint coordinate system for the clinical description of three-dimensional motions: application to the knee. J Biomech Eng 1983;105:136–44. [30] Harrison MM, Cooke TD, Fisher SB, Griffin MP. Patterns of knee arthrosis and patellar subluxation. Clin Orthop Relat Res 1994;5:6–63. [31] Spoor CW, van Leeuwen JL. Knee muscle moment arms from MRI and from tendon travel. J Biomech 1992;25:201–6. [32] Kim HJ, Fernandez JW, Akbarshahi M, Walter JP, Fregly BJ, Pandy MG. Evaluation of predicted knee-joint muscle forces during gait using an instrumented knee implant. J Orthop Res 2009;27:1326–31. [33] Andriacchi TP, Mundermann A, Smith RL, Alexander EJ, Dyrby CO, Koo S. A framework for the in vivo pathomechanics of osteoarthritis at the knee. Ann Biomed Eng 2004;32:447–57. [34] Goudakos IG, Konig C, Schottle PB, Taylor WR, Singh NB, Roberts I, et al. Stair climbing results in more challenging patellofemoral contact mechanics and kinematics than walking at early knee flexion under physiological-like quadriceps loading. J Biomech 2009;42:2590–6. [35] Heegaard J, Leyvraz PF, Van Kampen A, Rakotomanana L, Rubin PJ, Blankevoort L. Influence of soft structures on patellar three-dimensional tracking. Clin Orthop Relat Res 1994;23:5–43. [36] Cohen ZA, McCarthy DM, Kwak SD, Legrand P, Fogarasi F, Ciaccio EJ, et al. Knee cartilage topography, thickness, and contact areas from MRI: invitro calibration and in-vivo measurements. Osteoarthritis Cartilage 1999;7: 95–109. [37] Carter DR, Wong M. Modelling cartilage mechanobiology. Philos Trans R Soc Lond B Biol Sci 2003;358:1461–71. [38] Gorniak GC. Patterns of patellofemoral articular cartilage wear in cadavers. J Orthop Sports Phys Ther 2009;39:675–83.