6 Material and methods
6.1 Protocols, model representations, material properties, and
Study I: Knee joint imaging and segmentation: A previously described numerical model was applied and modified in this analysis [216]. An asymptomatic male subject (28 years, 80 kg) was imaged in an MRI scanner and the knee joint tissues (femoral, tibial and patellar cartilages, and meniscus) were segmented and meshed. The anatomical geometries were segmented using the MRI data (a clinical 3.0 T scanner, Philips, Best, Netherlands) which were taken using a 3D FSE sequence (VISTA) (in-plane resolution = 0.5 mm, slice thickness = 0.5 mm, TR = 1300 ms, TE = 32.3 ms). This MRI acquisition was conducted according to the ethical guidelines of Kuopio University Hospital, Finland. This imaging procedure was conducted with the permission (94/2011) from the local ethical committee of the Kuopio University Hospital, Kuopio,
Study Geometry Study aim
Constitutive material Cartilage Menisci Ligaments
I 3D knee joint
Effect of different constitutive representations and structural constituents of ligaments on cartilage. cartilage degradation model in injured tissue.
Finland. In addition, patellar and quadriceps tendons (PT and QT), as well as six knee ligaments (ACL, PCL, MCL, LCL and medial and lateral patellofemoral ligaments) were segmented and meshed in Abaqus v.6.14 (Dassault Systémes, Providence, USA).
Knee tissue representations and biomechanical properties: To evaluate the effect of ligament geometry and the contributions of their structural constituents (fibril network, nonfibrillar matrix, fluid) to the knee joint response during the stance phase of gait, springs and 3D continuum geometries were implemented in FE knee models.
In general, cartilages and menisci were meshed with hexahedral pore pressure elements (type = C3D8P). In the simplified case, the ligaments and tendons were modeled using an array of spring elements [216,253]. In the more sophisticated case, 3D continuum geometries for ligaments were meshed with tetrahedral pore pressure elements (type = C3D4P). For the spring ligaments, a pre-tension at their segmented length in MRI data was selected (5% on ACL and PCL and 4% on LCL and MCL [215,216,254]) while in the 3D continuum ligaments, pre-strains were imposed utilizing the predefined field stress (the equivalent stress tensor to the pre-strain was calculated and applied for each ligament). In total, five knee joint models were developed with different constitutive material formulations for the ligaments: i) spring, ii) linear elastic, iii) hyperelastic, iv) porohyperelastic and v) fibril-reinforced porohyperelastic (FRPHE) material. A group of material constants for each constitutive model was estimated for each ligament based on experimental studies (ACL [3,209,219,255–258,258–263], PCL [3,75,209,255,257,260,262,264–266], MCL [3,77,106,210,218,219,267–272], LCL [3,81,106,217,268,270,271,273–278], PT [261,268,273,279–291], and QT [261,289,290,292]) (Table 6.2).
Figure 6.1: General overview of the experimental and computational in vitro and in vivo methods used in the thesis.
Table 6.2: Reference values for ligament material properties (study I).
k = spring constant, E =Young's modulus, ν = Poisson's ratio, C1 = material coefficient, D = compressibility coefficient, k0 =initial permeability, M = exponential term for the strain-dependent permeability, Em = nonfibrillar matrix modulus, Ef = fibril network modulus, νm = Poisson’s ratio of the nonfibrillar matrix, C = density ratio between primary and secondary fibrils. †Adjustment method is explained in section 6.2.
ACL [3,209,219,255–258,258–263], PCL [3,75,209,255,257,260,262,264–266], MCL [3,77,106,210,218,219,267–
272], LCL [3,81,106,217,268,270,271,273–278], PT [261,268,273,279–291], and QT [261,289,290,292]).
On the other hand, for all FE knee models, cartilages were defined as an FRPVE material [58,195,241], while menisci were considered as an FRPHE model. Material properties have been validated and tested in previous FE modeling studies [253,293]
(Table 6.3).
Boundary conditions: The implementation of the subject´s gait information and the determination of boundary conditions in the FE knee model were conducted using a similar approach as in previous studies [21,159,215,216,294,295]. Briefly, the translational forces (inferior–superior, medial-lateral and anterior-posterior), moments (in the internal-external and valgus–varus directions) and extension-flexion rotation of the subject were applied to the reference point, which is located at the center of the transepicondylar axis of the femur. In addition, the quadriceps force was implemented in the QT and followed the direction of the distal-proximal axis of the femur, and distal nodes of the tibia were defined fixed during the motion. The cartilage-bone interfaces were considered as rigid. Frictionless surface-to-surface contact was established for all the surfaces in contact, including the contact interactions between cartilages and menisci as well as the external surfaces of ligaments and cartilages, and PCL and ACL (Figure 6.2).
Table 6.3: Material parameters implemented for cartilages and menisci (study I).
𝐸0 = initial fibril network modulus, 𝐸f = fibril network modulus, 𝐸ε = strain-dependent fibril network modulus, η = damping coefficient, 𝐸nf = nonfibrillar matrix modulus, 𝑘0 = initial permeability, 𝜈nf = Poisson´s ratio, M = exponential term for the strain-dependent permeability, C = ratio of primary to secondary collagen fibers, nfl = depth-wise fluid fraction distribution, h indicates normalized distance from the cartilage surface (surface = 0, bottom = 1).
*Julkunen et al. [296]
**Wilson et al. [195]
†Saarakkala and Julkunen [297]
††Dabiri and Li [298]
‡Makris et al. [68]
Material
parameter Femoral cartilage Tibial cartilage Patellar cartilage Menisci
𝐸0 (MPa) 0.92* 0.18* 1.88* -
𝐸f (MPa) - - - 28††
𝐸ε (MPa) 150* 23.6* 597* -
η (MPa) 1062* 1062* 1062* -
𝐸nf (MPa) 0.215* 0.106* 0.505* 0.5*
𝑘0 (10-15 m4N-1s-1) 6* 18* 1.9* 1.25††
𝜈nf (-) 0.15** 0.15** 0.15** 0.36††
M (-) 5.09* 15.64* 15.93* 5.09
C (-) 12.16** 12.16** 12.16** 12.16**
nfl (-) 0.8 – 0.15h** 0.8 – 0.15h** 0.8 – 0.15h** 0.72‡
1
Figure 6.2: General overview of study I. Workflow used for generation of tissue geometries and implementation of subject´s gait data into the FE model.
Study II: Injurious compression experiments: Explants of articular cartilage were harvested from the patellofemoral grooves of 1-2 week-old calves. Full-thickness tissue cylinders were obtained using a 3-mm dermal punch, and the top 1-mm disk was harvested with a razor blade [299]. Disks of articular cartilage were incubated in serum-free medium and equilibrated for 2 days before the experiments. Treatment groups were matched for the location of the cartilage cylinders along the knee joint [125,300]. On day 0, injurious compression experiments were performed on twenty two cartilage explants using a custom-designed, incubator-housed loading device [125,301]. Each cartilage explant was subjected to unconfined compression at a strain rate of 100%/s to 50% final strain, followed by a release at the same rate. Then, twelve disks from those twenty two samples were analyzed at day 0 and the remaining explants (n = 10) went to the dynamic loading protocol (Figure 6.3).
Cartilage cylinders were placed in the loading chamber and then subjected to dynamic compression after the injurious compression on day 0. Cartilage disks were subjected to unconfined dynamic compression (haversine waveform with 15%
amplitude [125]) using a displacement-controlled loading (1 Hz) continuously for 1 h, followed by a 5 h rest period with the applied load removed. This loading cycle was repeated 4 times per day for 12 days. The control group (no injury + 15% strain) was subjected to the same protocol, however there was no initial injurious loading.
Thus, a dynamic compression at 15% strain was conducted on two different groups of explants at the same time: i) 10 explants from injurious loading followed by 15%
dynamic strain, ii) 11 explants from no injurious loading followed by 15% dynamic strain. In addition, 19 extra samples were immersed in the same culture medium, but no dynamic loading was applied with 10 of those disks being studied on day 0 (no injury group) and nine samples on day 12 (free swelling group) (Figure 6.3).
The mechanical parameters for the numerical models were determined based on three additional intact cartilage explants placed in the loading chamber and compressed in the incubator-housed loading device [299]. Initially, a pre-strain of 10% of the cartilage thickness was applied (200 s compression, followed by 600 s relaxation). Then two additional stress-relaxation steps of 2.5% strain, each, were imposed (30 s compression, 300 s relaxation). The biomechanical parameters (𝐸f , 𝐸nf, and 𝑘) were obtained by fitting the numerical models in Abaqus to the unconfined compression experiments using Matlab by minimizing the normalized mean squared error between experimental and numerical forces [302] (fminsearch function) (Table 6.4).
Digital densitometry and cell viability analysis: After the dynamic loading measurements, histological sections of the samples were then cut perpendicularly to the cartilage surface. The depth-wise optical density profiles were quantified utilizing a digital densitometry (DD) method [303] for the five groups. Microscopic sections were stained with Safranin-O under standardized conditions (pH 4.6). The Safranin-O stain amount in the histological sections was determined by capturing grayscale images from the sections with a light microscope (Olympus CH-2, Olympus, Tokyo, Japan). Three optical density image maps from three independent sections of each sample were examined and averaged perpendicular to the depth-wise direction to compute a depth-depth-wise optical density profile. Finally, the profiles of the depth-wise optical density were utilized for estimating the depth-wise distribution of the FCD content [304,297,305].
Cell viability experiments were performed on day 0 and 12 days following the injury to evaluate the effect of dynamic loading on local chondrocyte viability. One mm-thick slices were collected from 1-2 additional samples in each group and then assayed and imaged with two fluorescence microscope filters to generate an image which was used in the viability analysis. Viable green cells were stained using
fluorescein diacetate while propidium iodide (both from Sigma-Aldrich, St. Louis, MO, USA) stained dead cells red [125,306].
Figure 6.3: General overview of study II. Schematic representation of the custom-designed apparatus and timeline used to perform injurious dynamic compression experiments.
Measurements on day 0 and day 12 were made for separate groups of cartilage explants.
Schematic representation of the experimental design and number of samples in each group.
Groups with stars and polygons indicate histology and cell viability analysis, respectively.
Representative histological images were used to obtain the optical density distribution and create FE element geometries.
Numerical modeling: In the mechanobiological model, two-dimensional (2D) FE models were developed for three representative samples with cracks in the cartilage examined after injurious compression. Images from histological samples were imported into Mimics v19.1 (Materialise, Belgium) where segmentation of the tissue lesions was done. Then, the segmented geometries were imported and meshed with pore pressure continuum elements (type CPE4P) in Abaqus. A FRPHES constitutive model [241] was defined to model the articular cartilage samples. The depth-dependent collagen density content and orientation of the primary fibrils were varied and then implemented into the models based on experimental studies of bovine calves [299,307,308]. The initial depth-wise profiles of the FCD were obtained as an average of the estimated optical density profiles of 12 cartilage samples. The depth-wise optical density profile was transformed into the FCD content profile by using the mean value of the FCD (in mEq/ml) reported for calf femoral cartilage [57].
Table 6.4: Material parameters for the FRPES cartilage model (study II).
*Obtained by fitting the model to unconfined compression experiments.
**Obtained from Li et al. [309], the effective Poisson’s ratio of the tissue became close to 0.1 due to reinforcement by the collagen fibrils.
***Obtained from Wilson et al. [241]
†Obtained from optical density distributions and converted into FCD values by using the average value of FCD reported by Han et al. [57]
††Obtained from Saarakkala and Julkunen [297]
h indicates normalized distance from the cartilage surface (surface = 0, bottom = 1).
FRPES = fibril-reinforced poroelastic swelling.
Boundary conditions: Following an initial free-swelling step, the FE explant model was subjected to the unconfined dynamic compression protocol as described in the experimental section. Briefly, the bottom of the sample was fixed in the axial direction and fluid flow was allowed through free and lesion surfaces [21,25] while the contact
Material
parameter Value Description
𝐸f (MPa) 20* Fibril network modulus
𝐸nf (MPa) 0.16* Nonfibrillar matrix modulus
𝑘 (10-15 m4N-1s-1) 1.3* Hydraulic permeability
between the plate and cartilage explant was assumed as being impermeable and frictionless.
Cartilage degeneration algorithm: To simulate the degradation of injured cartilage under relevant mechanical loading, an iterative degeneration algorithm was developed which is associated with three different outputs: interstitial fluid velocity, maximum shear, and deviatoric strains. These outputs were transferred to the algorithm via Abaqus-Matlab. The algorithm assumes that an FCD loss can occur if either the deviatoric strain (𝜀dev) is greater than a threshold of 20%, or the maximum shear strain (𝜀shr) exceeds 50% or fluid velocity (𝑣fl) is more than 0.04 mm/s during the complete loading cycle. The fluid velocity, as well as deviatoric and maximum shear strains, were defined previously by Equations 4.13, 4.28, and 4.29, respectively.
In particular, the criteria related to the deformation were defined based on previous investigations [22,136,310,311] and an additional sensitivity analysis (for both, shear and deviatoric strain threshold range was set from 0.15 to 0.7) while the criterion regarding fluid velocity threshold (0.04 mm/s) was established with reference to the preliminary analysis where the initial fluid velocity value was obtained from the intact cartilage model (Figure 6.4).
In addition, different linear and nonlinear approaches to simulate the time-dependent degradation evolution were tested; these showed similar results. Thus, a piece-wise constant degeneration rate factor 𝐷r was implemented which can be defined as
𝐷r= {
0 if 𝜀𝛽< 𝜀𝛽,thres,
0.6 if 𝜀𝛽,thres≤ 𝜀𝛽 ≤ 𝜀𝛽,failure , 𝛽 = dev, shr 1 if 𝜀𝛽> 𝜀𝛽,failure,
(6.1)
where 𝜀𝛽 is the strain variable (either deviatoric (β = dev) or maximum shear strain (β = shr)), 𝜀𝛽,thres is the strain threshold at which the non-fibrillar matrix damage is initiated and 𝜀𝛽,failure is the value for a complete failure of the ground substance (𝜀𝛽,failure = 1). Likewise, a similar approach was used in the case of the fluid velocity driven degeneration
𝐷r= {
0 if 𝑣fl< 𝑣fl,thres, 0.6 if 𝑣fl,thres≤ 𝑣fl≤ 𝑣fl,failure , 1 if 𝑣fl> 𝑣fl,failure,
(6.2)
where 𝑣fl,thres and 𝑣fl,failure are the fluid velocity values at which the non-fibrillar matrix damage is initiated (𝑣fl,thres = 0.04 mm/s) and the entire failure occurs (𝑣fl,failure
= 1.5 mm/s), respectively.
Non-fibrillar matrix degradation (PG depletion) was modeled based on the FCD decrease method, described previously by Equation 4.32, considering the effect of the degeneration rate factor 𝐷r. After every iteration, a new FCD content was implemented and subsequently, a new loading cycle was simulated (Figure 6.4). In order to avoid computational instabilities, the FCD content was not allowed to reach zero during the degenerative iterations. Thus, the minimum allowed FCD content was set to 10% of the smallest initial FCD content, which was obtained from the superficial layer.
Figure 6.4: General workflow of the iterative algorithm used to evaluate the degeneration mechanisms. FCD initial content, sample geometries, and mechanical properties were obtained from experimental measurements; collagen arrangement, depth-wise collagen distribution and osmotic properties were taken from the literature. After the mechanical simulation, fluid velocity, maximum shear strain and deviatoric strains were obtained and utilized for driving the cartilage softening algorithm, and new FCD content was implemented into the FE model. This loop was repeated 50 times, and the FCD predictions were compared with the FCD decrease obtained from DD experiments. As convergence criterion, the FCD content distribution was ensured to reach an equilibrium state, showing null FCD loss progression.
Study III: MRI imaging protocol: Two patients with ACL reconstructed knees (Patient 1: 44 years, 79 kg; Patient 2: 39 years, 62 kg) were imaged at different follow-up time points using an MRI scanner (3.0T, MR750w, General Electric Healthcare). The patients gave informed consent, and the obtained data were permitted by and carried out in agreement with the regulations of the Institutional Review Board (#11-06734) under the Human Research Protection Program at UCSF. The knee joint geometries were obtained with a 3D intermediate-weighted, fluid sensitive, fat-saturated fast-spin-echo (CUBE) [312]. The acquisition parameters were: TR = 1500 ms, TE = 25 ms,
field of view (FOV) = 16 cm, slice thickness = 1 mm, matrix = 384x384, pixel size = 0.42 mm x 0.42 mm. In addition, acquisition of T1ρ/T2 [31,313] was performed to quantify the T1ρ and T2 relaxation times in the knee joint at 1 and 3 years after surgery, using the following acquisition parameters: TR/TE = 9ms/min full, FOV = 14cm, matrix = 256x128, pixel size = 0.55 mm x 1.10 mm (0.55 mm x 0.55 mm with interpolation), slice thickness = 4mm, views per segment = 64, lock frequency = 500 Hz, time of spin-lock = 0/10/40/80ms for T1ρ and preparation TE = 0/13.7/27.3/53.7ms for T2 [314].
Furthermore, the T1ρ and T2 maps were computed on a pixel-by-pixel basis using a mono-exponential fit in Aedes (http://aedes.uef.fi) and in-house plugins for Matlab.
Patients’ gait data and musculoskeletal modeling: The gait information of the patients was obtained at the 1-year follow-up time point using a published protocol [315]. The motion capture system consisted of 10 cameras (Vicon, Oxford metrics) and 41 retro-reflective markers to obtain segment position data. Two in-ground force plates (AMTI, Watertown, USA) were utilized to obtain the ground reaction forces simultaneously. Then, to calculate knee joint contact and lower extremity muscle forces for the FE knee joint model, a generic musculoskeletal model in OpenSim (SimTK, Stanford, CA, USA) was used [316]. Anatomical landmarks collected during gait were used for scaling the generic model to the individual anthropometrics [172].
Knee joint and quadriceps forces were simulated for every trial. Then, similar to the methodology presented in study I, the average of the flexion-extension angles, joint moments, knee translational forces, and quadriceps forces during the stance phase were computed and implemented to drive the FE knee joint models (Figure 6.5).
FE knee joint geometries: The 3D CUBE MRI data obtained at the 1-year time point was used to segment the knee geometries (femoral, tibial, patellar cartilages, menisci, collateral, and cruciate ligament insertions), including the lesion/defect in the lateral tibial cartilage (Patient 1: depth = 3.3 mm, width =~ 2.8 mm) and in the lateral femoral cartilage (Patient 2: depth = 1.1 mm, width =~ 1.6 mm). Then, the geometries were imported and meshed in Abaqus where the FE models were developed. Cartilages and menisci were meshed using 8-node hexahedral poroelastic (C3D8P) elements and modeled using FRPVES and FRPES materials, respectively. Material parameters used in study III are given in Table 6.5.
Boundary conditions: Similar approach used in study I, the tibial cartilage-bone interface was fixed in all directions and fluid flow was allowed only through the inner lesion surfaces (zero fluid pore). After an initial free-swelling step, the stance phase of the patient’s gait obtained from the musculoskeletal model was implemented identically with earlier studies [216,294,295,316].
Figure 6.5: General overview of study III. Knee geometries were generated and acquired from the T1ρ and T2 maps. Gait data was used in the musculoskeletal model to obtain the joint moments, extension-flexion angle and knee contact forces. FE joint models were subjected to dynamic loading (stance phase of gait). After the simulation, maximum shear strain, deviatoric strain, and fluid velocity were obtained and used for driving the cartilage degenerative algorithm, and new fixed charge density (FCD) content was implemented in the mechanobiological finite element model. FCD predictions were contrasted with the T1ρ and T2 maps for verifying the progression of post-traumatic OA.
Table 6.5: Material parameters implemented for cartilages and menisci (study III).
𝐸0 = initial fibril network modulus, 𝐸f = fibril network modulus, 𝐸ε = strain-dependent fibril network modulus, η = damping coefficient, 𝐸nf = nonfibrillar matrix modulus, 𝑘0 = initial permeability, 𝜈nf = Poisson´s ratio, M = exponential term for the strain-dependent permeability, C = ratio of primary to secondary collagen fibers, nfl = depth-wise fluid fraction distribution, cFCD = depth-wise fixed charge density distribution, ρz = depth-wise collagen distribution, h indicates normalized distance from the cartilage surface (surface = 0, bottom = 1). The first set of mechanical properties are identical to those of study I.
*Julkunen et al. [296]
**Wilson et al. [195]
†Saarakkala and Julkunen [297]
††Dabiri and Li [298]
‡Makris et al. [68]
‡‡Mow and Ratcliffe [47]
Mechanobiological knee joint models: The cartilage degeneration algorithm developed and presented in study II was implemented into the FE knee joint models. Iterative simulations using this validated degradation algorithm combined with the 3D FE knee models were carried out to predict the variations in FCD content in injured cartilage. As described previously, this algorithm has been developed to be driven by three different model outputs: deviatoric and maximum shear strains, and interstitial fluid velocity. Since the differences in the inherent tissue mechanical properties regarding those analyzed in study II, the mechanobiological model assumed FCD loss to occur in the aforementioned FRPVES material if either the deviatoric strain (𝜀dev) exceeds a threshold of 20%, or the maximum shear strain (𝜀shr) is greater than 40%, or the fluid velocity (𝑣fl) is larger than 0.03 mm/s during the entire stance phase of the gait cycle. These thresholds were selected based on further sensitivity analysis conducted to determine the link between numerical results and the changes observed in the MRI maps.
Material
parameter Femoral cartilage Tibial cartilage Patellar cartilage Menisci
𝐸0 (MPa) 0.92* 0.18* 1.88* -
Comparison between FCD decrease predictions and MRI maps. The FCD loss predictions from mechanobiological models were contrasted with the changes observed in T1ρ and T2 relaxation times of cartilage from 1 to 3 years after surgery in both qualitative and quantitative manner. In the quantitative analysis, the degenerated volumes were calculated from patient-specific FE models and MRI follow-up maps. The degenerated volume from MRI was defined as the volume of cartilage with either T1ρ
Comparison between FCD decrease predictions and MRI maps. The FCD loss predictions from mechanobiological models were contrasted with the changes observed in T1ρ and T2 relaxation times of cartilage from 1 to 3 years after surgery in both qualitative and quantitative manner. In the quantitative analysis, the degenerated volumes were calculated from patient-specific FE models and MRI follow-up maps. The degenerated volume from MRI was defined as the volume of cartilage with either T1ρ