Studies I-IVuse MRI and motion capture collected at the University of California San Francisco (UCSF). All subjects included in this study gave informed consent and data acquisition was approved by and carried out in accordance with the rules and regulations of the Institutional Review Board under Human Research Protection Program at UCSF. An overview of the subject demographics is provided inTable 6.1. To generate subject-specific FE models, the following inputs were required: knee joint geometry, knee joint motion, meshing and selection of material formulation for the soft tissues.
Table 6.1: Subject characteristics at 1-year follow-up time-point inStudies I-IV.
Study I Study II Study III Study IV
Group ACLR ACLR ACLR Control ACLR
Subjects (n) 1 2 7 6 2
Age (years(std)) 34(-) 41.5 (3.53) 37.8 (6.5) 31.3(0.37) 41.5 (3.53)
Gender (m/f) 0/1 2/0 5/2 3/3 2/0
Weight (kg(std)) 63.6(-) 70.5(12.0) 70.6(13.4) 68.3(10.3) 70.5(12.0) Height (m(std)) 1.6(-) 1.8(0.02) 1.74(0.1) 1.69(0.1) 1.8(0.02)
Knee joint geometry
In all studies 3D fast spin-echo sequence (FSE) MR images were acquired at the 1-year post-ACLR time-point using a 3T MR scanner (MR750w, Widebore) and a 1Tx/8Rx knee coil at UCSF. The 3D-FSE sequence is an intermediate-weighted, fluid-sensitive, fat-saturated sequence (CUBE sequence). The acquisition parameters were: repetition time (TR) = 1500ms, echo time (TE) = 25ms, field-of-view (FOV) = 16 cm, slice thickness = 1 mm, echo train length (ETL) = 32, matrix = 384x384 and pixel size = 0.42×0.42 mm2. Then, Seg3D (v2.4.4,
Figure 6.2: An example of the segmentation and post-processing steps. a) Output from manual segmentation software; b) Corresponding mesh for surface in a); c) Smoothing and reducing mesh density step; d) Final mesh used in FE models.
NIH/NIGMS CIBC, USA) was used to manually segment the knee joint soft tissues from the 3D CUBE MR images at 1-year after ACLR (Fig. 6.2a). After post-processing in MeshLab (v2016.12 [263]) and MATLAB (v2016b and 2018b, MathWorks, USA) (Fig. 6.2b and c), the solid geometries were imported into
Abaqus (v6.13-3 and v2018, Dassault Systemes, USA) for final meshing (Fig. 6.2d) and model creation (Fig. 6.3). Two geometries (in terms of model complexity) were created from the manually segmented MR images. The complex geometry (Geometry A, Fig. 6.3) included tibial, femoral and patellar cartilage, menisci, collateral and cruciate ligaments, patello-femoral ligaments, quadriceps and patellar tendons. Geometry A was used in Studies I and IV. The simplified geometry (Geometry B,Fig. 6.3) only contained tibial cartilage, femoral cartilage, menisci, cruciate and collateral ligaments. Geometry B was used inStudies I-III.
Two patients with ACLR had focal full-depth cartilage defects (Fig. 6.4). These
Figure 6.3: Overview of the geometry generation and schematic presentation of the differences between complex geometry (Geometry A) and simple geometry (Geometry B) utilized in the current thesis and created based on manual segmentation.
patients were used inStudies II-IV. These lesions were manually segmented and included in the FE analysis. For Patient 1, the defect was located on the lateral tibial cartilage, with a depth of 3.3 mm and a width of 2.8 mm. ForPatient 2, the lesion was on the lateral femoral cartilage, with a depth of 1.1 mm and a width of 1.6 mm. To ensure the accuracy of the segmentation process, we consulted an experienced musculoskeletal radiologist.
Figure 6.4: MRI and segmented 3D geometry of the focal cartilage defects forPatient 1(upper row) andPatient 2(lower row).
Geometry mesh
In all studies, 8-node hexahedral poroelastic (C3D8P) elements were used for articular cartilage and corresponding elastic elements (C3D8) were used for the menisci. The characteristic elements length of femur, tibial and patellar cartilage were 1.5, 0.5, 0.7 mm respectively. For menisci, the characteristic element length was 1 mm. For all soft tissues three depth-wise element layers were enforced. The characteristic element lengths were obtained from previous studies [230, 234], where mesh convergence tests were conducted to ensure that the finite element analysis results would not be influenced by changes in mesh density (i.e. changes in element size) for the tibial cartilage stresses, strains and pore pressures.
Knee joint motion
Subject gait data was measured at UCSF using a previously established protocol [180]. For each subject, position data was recorded at 250 Hz using a 10-camera motion capture system (Vicon, UK). The ground reaction forces were measured with two in-ground force plates (AMTI, USA) at 1000 Hz. Segment position data was collected using 41 retro-reflective markers. For each subject, four gait trials were measured at a controlled walking speed of 1.35 m·s−1. Afterwards, a rigid-body seven-segment subject-specific musculoskeletal model was developed in Visual3D (C-Motion, USA). Joint rotations were calculated using Cardan sequence, while joint moments were calculated using standard inverse dynamics.
More details about the marker positioning and the musculoskeletal model are provided in Samaan et al. [180]. In all studies the knee joint motion, obtained from motion capture data, was taken as the average of three gait trials (Fig. 6.5). In Studies I-III, the knee joint translational forces and the knee joint rotations were directly used as inputs for the FE models. In Studies I and IV, the knee translational forces and the joint moments were used as inputs for the FE models.
InStudies I and IV, the quadriceps muscle force, divided into anterior-posterior
Figure 6.5: Translational forces, knee joint rotations and knee joint moments obtained from motion capture.
and distal-proximal components, was computed in OpenSim (v. 4.0, SimTK, USA).
The generic Gait2392 [238] model was scaled in order to match the generic marker
positions with the marker positions from Visual3D. The inverse kinematics data from Visual3D was used as input for the musculoskeletal simulation in OpenSim.
Then Residual Reduction Algorithm (RRA), a form of forward dynamics, was coupled with Computed Muscle Control (CMC). RRA minimizes large non-physical compensatory forces (i.e. residuals) by altering the torso mass center in the model, thus ensuring dynamic consistency with ground reaction force data [264, 265]. CMC is used to compute lower extremity muscle forces using static optimization, by reducing muscle redundancies and includes muscle activation and contraction dynamics [265, 266]. The OpenSim-CMC approach has been shown to produce similar muscle activation patterns as electromyography measurements (EMG) [242, 267, 268]. The quadriceps muscle force was then calculated as the algebraic sum of the rectus femoris, vastus lateralis, intermedius and medialis muscles [205].
FE models
Figure 6.6: Overview of geometries and motion inputs for modelsA,B,CandD.
Four FE models were created in this thesis (Fig. 6.6):
• Model A uses kinetic actuation and Geometry A. This model is the most complex used in this thesis and similar to the one developed by Halonen et al. [205]. The implemented moments included internal-external and varus-valgus, with flexion-extension as rotation. The moment scaling was 50% of the measured values, similarly as before [205], indicating that muscles absorb half of the measured moments. By this assumption, a good match between modeled and experimental femoral rotations was obtained.
Translational forces (anterior-posterior, medial-lateral and distal-proximal) and quadriceps muscle forces (divided in anterior-posterior and distal-proximal components) were implemented in this model. We refer to ModelAalso as the kinetic driven model, even though the flexion-extension angle was implemented with kinematics.
• ModelBuses kinetic-kinematic actuation andGeometry A. The implemented translational and quadriceps forces were the same as in Model A, while flexion-extension and internal-external motions were implemented as rotations instead of moments. The varus-valgus rotation was kept free, since the variation in motion capture for varus-valgus rotation is known to be higher than with the other two rotation directions [235–237]. This assumption has been made in several earlier studies [205, 230].
• ModelCuses kinetic-kinematic actuation andGeometry B. For ModelC, only the translational forces and rotations were implemented.
• ModelDuses kinetic-kinematic actuation andGeometry B. In addition to the translational forces and rotations, the anterior-posterior component of the quadriceps force was added to the anterior-posterior component of the translational forces. This was done in order to consider the quadriceps force contribution without a need for more complex geometry with patella and quadriceps and patellar tendons.
The models were used as follows:
• Models A-D were used in Study I (n=1) to compare different motion implementation methods.
• ModelCwas used in StudyII(n=2) andStudy III(n=13) to rapidly assess the risk of OA onset and development due to collagen degeneration and/or PG loss.
• Model Awas used in conjunction with an adaptive degeneration algorithm in Study IV (n=2) to simulate the changes in FCD content around cartilage lesions.
6.2.2 Material properties
Several studies have shown that simpler material formulations for articular cartilage [190], menisci and ligaments [204] can produce similar, if not the same response as more complex material formulation. These simpler formulations are also less computationally intensive and may be better suited for clinical applications. In Studies I-III, articular cartilage was modeled as TIPE [190], and
menisci were modeled as TIE [33, 67, 190, 191]. The material properties for cartilage and menisci are shown inTable 6.2.
Ligaments were assumed to be in tension at their normal length determined from MRI, and therefore a prestrain of 5% was applied to ACL and PCL, while 4%
was applied to MCL and LCL [269–271]. Similar as before [205], the stiffness values were: k=120 N/mm for healthy ACL; k=200 N/mm for PCL; k=100 N/mm for MCL and LCL. For patients with ACLR, the graft used in the reconstruction was a posterior tibialis tendon with a stiffness 380 N/mm [271]. The meniscal attachments were modeled as linear springs with a total spring stiffness of 350 N/mm for each meniscal horn [272]. For Models A and B, the stiffness of the quadriceps and patellar tendons were 475 N/mm and 545 N/mm, respectively [273]. Additionally, the lateral and medial patellofemoral ligaments (LPFL and MPFL) were modeled as linear truss elements [199, 270]. The elastic modulus was defined as 17 MPa for LPFL and 19 MPa for MPFL and the Poisson’s ratio was set to 0.499 for both [205, 274].
Table 6.2: Material parameters of cartilage, meniscus used in the FE models in Studies I-III.
TIPE / TIE Ep Et νp νtp Gt k(10−15 e0
(MPa) (MPa) (−) (−) (MPa) m4/Ns) (−)
Cartilage 24 0.46 0.42 0.06 12 1 4
Meniscus 20 159.6 0.3 0.01 50 -
-Parameters: Ep– in-plane Young’s modulus, Et – out-of-plane Young’s modulus, νp– in-plane Poisson’s ratio,νtp– out-of-plane Poisson’s ratio,Gp– in-plane shear modulus,Gt– out-of-plane shear modulus, k– permeability,e0– initial void ratio.
Note:Gtis calculated directly in Abaqus usingGp=Gp/2(1+νp). ACLR ligament parameters were used for patients with reconstruction; ACL ligament parameters were used for controls.
In a previous study by Orozco et al. [203], a novel mechanobiological FE model was developed to predict FCD loss in mechanically injured cartilage. The iterative degeneration algorithms proposed in that study were also used in Study IV.
Similar to Orozco et al. [203], articular cartilage was modeled as fibril reinforced poroviscoelastic (FRPVE) with swelling and menisci as fibril reinforced poroelastic (FRPE) materials [275]. Ligaments and tendons were modeled identically as in Studies I-III. The depth-dependent primary collagen fibril orientations with split-line patterns, fluid fractions and FCD distribution were implemented in the cartilage tissue. For menisci, the primary collagen fibrils were oriented circumferentially, and the fluid fraction and FCD content were homogeneously distributed. The methodology is based on previous studies [204, 276]. A summary of the material properties is shown inTable 6.3.
6.2.3 Boundary conditions and simulation steps
Boundary conditions: The FE models used three reference points (Fig. 6.3): RP -femoral reference point, taken as the midpoint of the distance between medial and lateral femoral epicondyles [230, 278]; RPPatella - patellar reference point, taken as
Table 6.3: Material parameters of cartilage, meniscus and ligaments used in the FE models ofStudy IV.
Material Cartilage Menisci
parameter Femur Tibia Patella
E0(MPa) 0.92a 0.18a 1.88a
-Ef(MPa) - - - 28b
Ee(MPa) 150a 23.6a 597a
-η(−) 1062a 1062a 1062a
-En f(MPa) 0.215a 0.106a 0.505a 0.5a
k0 6a 18a 1.9a 1.25b
νn f(−) 0.15 [67] 0.15 0.15 0.36
M 5.09 15.64 15.93 5.09
C 12.16c 12.16c 12.16c 12.16c
Composition
nf l 0.8-0.15hc 0.72d
cFCD(mEq/ml) -0.1h2+0.24h+0.056e 0.03f
ρz 20.6h6-64.4h5+78+1h4-45.9h3+13.4h2-1.6h+0.96e 1
E0 - initial fibril network modulus; Ef - fibril network modulus; Ee - strain-dependent fibril network modulus; η - damping coefficient; En f - non-fibrillar matrix modulus;k0- initial permeability;νn f - non-fibrillar matrix Poisson’s ratio;
M- Exponential term for the strain-dependent permeability;C- ratio of primary to secondary collagen fibres; nf l - depth-wise fluid fraction;cFCD - depth-wise fixed charge density distribution; ρz - depth-wise collagen distribution;h - normalised distance from the cartilage surface (surface = 0, bottom = 1).
aJulkunen et al. [223],bDabiri et al. [258],cWilson et al. [67],dMakris et al. [78],
eSaarakala and Julkunen [56], fMow and Ratcliffe [277]
the centre of the patellar bone [234];RPQuad - quadriceps reference point, taken as the centre of the quadriceps tendon insertion point to the quadriceps muscles [234].
In all models (Models A - D), the knee joint translational forces, rotations or moments were implemented to the femoral reference point RP. Additionally, for the complex models (Models Aand B), the anterior-posterior and distal-proximal components of the quadriceps force were implemented into RPQuad). The cartilage-tibia bone interface was fixed for all six degrees-of-freedom (DOF) [230].
The meniscal and ligament insertion points in the tibial bone were also fixed for all six DOF [230]. Additionally, for Models A and B, the patellar tendon insertion point into the tibia was kept fixed [234]. The ligament insertion points into femur were kinematically coupled withRP. The tendon insertion points into the patellar bone were kinematically coupled with RPPatella [234]. The LPFL and MPFL insertion points into patella and femur were kinematically coupled with RP and RPPatella, respectively. Lastly, both RPand RPPatella (Master node) were tied to the respective cartilage surfaces (Slave surface) using∗T IE constraint in Abaqus. The contact between soft tissues was modelled using hard pressure-overclosure and frictionless surface-to-surface with node discretization [216, 230].
Simulation steps: Similar to previous studies [190, 204, 216, 230, 234], three steps were needed for the simulations:
1. For all models, the femoral cartilage was brought into light contact with the tibial cartilage and the menisci. In ModelsAandB, additionally the patellar cartilage was brought into contact with the femoral cartilage.
2. The femoral cartilage was rotated from its position in the MR image (i.e. ∼6◦ flexion angle, 0◦external angle) to the initial position at the beginning of the stance phase measured from motion capture (i.e. ∼2◦ flexion and external angle). ForStudies IIandIII, the varus-valgus rotation was also matched with motion capture information. Then, all initial forces (Models A-D), rotations (ModelsB-D) and moments (A) were implemented throughRP.
3. In all models, the time-dependent motion inputs were implemented.
For the analysis Abaqus/standard with soils consolidation (∗SOILS) dynamic analysis was used in all simulations.