6.4 Verification methods
6.4.2 T 1ρ and T 2 relaxation times
(6.2) whereεmax andεmin- maximum and minimum principal strains. Then a piece-wise constant degeneration rate factor (Dr) was defined as:
Dr =
(0 β<βThreshold
0.5 βThreshold≤β≤βFailure
1 β≥βFailure
(6.3)
whereβ,βthreshold andβf ailureare presented inTable 6.4.
Table 6.4: Thresholds and failure limits for deviatoric strains, absolute maximum shear strains and fluid velocity.
β βthreshold βf ailure
εdev(−) 0.2 1 εshr(−) 0.4 1 νf l(mm/s) 0.03 1.5
It should be noted that the deterioration rates provide predictions in a reasonable time and are in good agreement with experimental follow-up [203].
After each simulation, an updated FCD content was implemented into the FE model. Then a new computational iteration was executed, with the cartilage degeneration algorithm reiterated for 20 steps.
6.4 VERIFICATION METHODS 6.4.1 Literature
To investigate the FE model reliability inStudy I, the results were compared against experimental values reported in literature for knee joint rotations [92,282–287], knee joint reaction forces [16, 244, 288, 289] and contact pressures [216, 290].
6.4.2 T1ρ and T2 relaxation times
The T1ρ and T2 relaxation times were simultaneously acquired at both 1-year and 3-year follow-up time-points using a 3T MR scanner (GE Sigma HDx, General Electric Healthcare, USA) with a transmit/receive quadrature knee coil (Clinical MR Solution, USA) [26, 163] using a MAPSS sequence. The acquisition parameters were: TR/TE= 9 ms/min full, FOV= 14 cm, matrix= 256x128, pixel size= 0.55
mmx1.10 mm, slice thickness= 4 mm, views per segment= 64, spin-lock frequency=
500 Hz, time of spin lock (TSL)= 0/10/40/80ms for T1ρ and preparation TE=
0/13.7/27.3/53.7ms forT2. TheT1ρandT2relaxation times were then quantified by fitting the exponential decay functions (Eqs. 3.2 and 3.1, respectively) to the MRI signal. For theT1ρ relaxation times, the MRI signal acquired at the four TSLs was used. For theT2 relaxation times, the MRI signal acquired at the four preparation TEs was used. The T1ρ and T2 mapping was done using a two-parametric non-linear exponential fit with Aedes plugin for MATLAB and custom scripts, at both 1- and 3-year follow-up time points. Then the articular cartilage was segmented from theT1ρandT2mapped images using Aedes.
In Studies II-IV, the areas susceptible to OA onset and degeneration were compared against changes in T2 and T1ρ relaxation times both qualitatively and quantitatively. These are among the most established quantitative MRI parameters for articular cartilage, with T2 indicative of collagen network integrity and arrangement [154, 161, 162], and T1ρ of PG content [173, 174], respectively. In Studies III and IV, the results were also compared against WORMS grades evaluated at both 1- and 3-year follow-up time points from the high resolution MRI (please refer to the original publications in appendix).
Qualitative comparison - Axial views (Studies II-IV): Axial distribution maps were created for both the FE model parameters andT2 and T1ρ maps. For the FE models, the axial view shows the distribution of maximum principal stresses or absolute maximum shear strains on the articular cartilage surface (Fig. 6.7a). We used a custom Maltab script to determine theT2andT1ρvalues on the tibial cartilage surface (Fig. 6.7b).
Figure 6.7: Example of tibial surface definition for a) FE models and b)T2relaxation time maps.
Quantitative comparison -Predicted degeneration volumes (Studies II-IV): The total volumes susceptible to collagen and PG loss, predicted by the FE models at the 1-year follow-up time point, were compared against cartilage volumes with changes in the T2 and T1ρ relaxation times between 1- and 3-year follow-up time points.
Similarly, in StudyIV, the predicted volumes of FCD loss were compared against measured changes in theT2andT1ρrelaxation times, as follows:
1. For the maximum principal stress, absolute maximum shear strain and predicted FCD loss, we defined volumes-of-interest (VOI) for the medial and lateral tibial compartments. Where the respective thresholds were exceeded,
Figure 6.8: a) Overview of VOI (red) and total cartilage volume (blue) from FE model; b) Overview of VOI and total cartilage fromT2orT1ρmaps
the VOI was defined as the total volume of cartilage in which the thresholds were exceeded, expressed as a percentage of the total volume of the compartment (Fig. 6.8a).
2. For T2 and T1ρ maps, we defined VOIs for each compartment at both 1-year and 3-year follow-up times. The VOI was defined as the volume of articular cartilage with eitherT2orT1ρrelaxation time above 60 ms. Similar to step1, if no relaxation time exceeded this limit, the VOI was defined as “0”. Each VOI was expressed as a percentage of the total volume of the compartment both 1-year and 3-year follow-up time points (Fig. 6.8b).
3. We then subtracted the VOI at 3-year time point from the VOI at 1-year time point, reflecting the percentage of potentially damaged tissue from the total volume of each compartment.
Quantitative comparison -Correspondence between biomechanical parameters and quantitative MRI parameters(Study III): The relationship between biomechanical parameters (maximum principal stresses and absolute maximum shear strains) and changes in MRI parameters between 1- and 3-year follow-up time points was also quantified, as follows:
1. For FE models, the location and value of the peak maximum principal stress throughout the entire stance phase was measured for each joint compartment.
Similarly, the location and value of the peak absolute maximum shear strain throughout the entire stance phase was also measured for each joint compartment.
2. For MRI parameters, we evaluated the T2 relaxation time values at both 1-and 3-year follow-up times at the corresponding areas to those determined in Step 1 (defined as the average of 9 pixels). The change inT2relaxation time between the 1- and 3-year follow-up time points was computed. The change in T1ρ relaxation time between the 1- and 3-year follow-up time points was computed using the same procedure.
Statistical analysis (Study III): To quantify the relationship between the predicted and potentially damaged cartilage volumes we calculated Spearman’s correlation coefficients. Likewise, the Spearman’s correlation coefficients were calculated between the peak values of maximum principal stress and absolute
maximum shear strain and the change inT2and T1ρ relaxation times, respectively.
Bi-variate least square linear fits were also calculated.
There may be significant differences in joint mechanics between patients with ACLR and controls. This may lead to excessive tissue maximum principal stresses and/or shear strains in articular cartilage in patients with ACLR. Therefore, we evaluated differences in biomechanical parameters (peak values of maximum principal stresses and absolute maximum shear strains) and in MRI parameters (change inT2 and T1ρrelaxation times) between patients with ACLR and controls using Mann-Whitney U-test. We consulted an experienced statistician and all tests were carried out in MATLAB using custom scripts.
7 Results
7.1 METHODS TO IMPLEMENT MOTION IN FE MODELS
In Study I, all four models showed similar trends in varus-valgus rotation (Fig.
7.1a), all translation directions (Fig. 7.1c-e) and the medial-lateral joint reaction force (Fig. 7.1g). Differences between the kinetic and kinetic-kinematic driven models were found in the internal-external rotation (Fig. 7.1b), with an average difference over the stance phase of∼4◦. ModelsAandB, showed almost identical behaviour
Figure 7.1: Comparison between ModelsA, B,C and D in terms of tibiofemoral rotations, translations and reaction forces. The FE model results were compared against literature values for rotations [282–285, 291–294], translations [290–292] and distal-proximal reaction forces [16, 244, 288, 289].
in joint reaction forces through the stance phase. Model C showed the smallest
distal-proximal (Fig. 7.1f) and anterior-posterior (Fig. 7.1h) reaction forces, with the maximum differences of∼0.6 BW and∼0.25 BW, respectively, at∼20% of the stance when compared to ModelA. These differences reduced to∼0.05 BW for the rest of the stance phase. Even though ModelDonly included the anterior-posterior contribution of the quadriceps force, the maximum difference in the distal-proximal reaction force, as compared to ModelA, was reduced to∼0.25 BW (Fig. 7.1f). On the other hand, this force affected only slightly the anterior-posterior reaction force (Fig. 7.1h).
Importantly, the FE model results were within literature reported values for rotations [282–285, 291–294], translations [290–292] and tibiofemoral joint reaction forces [16, 244, 288, 289] (Fig. 7.1). Additionally, both the medial and lateral contact pressures were also within literature reported values [216, 290] (see original publication in Appendix).
In terms of maximum principal stresses (Fig. 7.2a, d), logarithmic strains (Fig. 7.2b, e) and pore pressures (Fig. 7.2c, f), in the medial compartment, the differences were up to two-fold between the kinetic driven model (Model A) and the kinetic-kinematic driven models (ModelsB,CandD) for all parameters at 10-30%
of the stance. Onward from∼30% of the stance, the differences in these parameters became very small. In the lateral compartment, all models behaved similarly, with the maximum differences of∼0.5 MPa,∼0.5% and∼1 MPa in maximum principal stresses, logarithmic strains and pore pressure, respectively, between ModelAand the rest of the models. The FE model results were also within literature reported
Figure 7.2: Comparison between Models A, B, C and D in terms of average maximum principal stress, average logarithmic strain and average pore pressure in the medial compartment (a, b and c) and lateral compartment (d, e and f).
values for the contact area [290, 295], contact pressure [233, 290], as well as,
maximum principal stresses, maximum principal logarithmic strains and pore pressure [213, 231, 233, 234] (seeFig. 7.3and the original publication in Appendix).
Figure 7.3: Medial and lateral contact areas (mm2) as a function of stance. Literature values were extracted from Gilbert2014 [290] and Zevenberg2018 [295].
For Models A and B, we also compared the patello-femoral joint kinematics and contact forces against literature reported values (Fig. 7.4). The patellofemoral translations and rotations were within literature reported values from a mobile x-ray imaging system [291]. The patellofemoral joint contact force was close to literature reported values for most of the gait cycle. [291, 296–298]
Figure 7.4: Comparison between ModelsAandBand literature values [291,297,298]
in terms of the patellofemoral kinematics (translations and rotations) and joint contact forces. Note: AP- anterior-posterior, IE – internal-external, ML – medial lateral, VV – varus-valgus, DP – distal-proximal, PF JCF – patellofemoral joint contact force.
Estimates of the total time needed for each model generation and simulation are shown inTable 7.1. In general, all models required almost the same run time.
However, the most complex, kinetic driven model Model (A) required several modifications in the contact controls between the femur and the patella. Thus, Model (A) needed 10 times longer to reach a converged solution, when compared to the simplest, kinetic-kinematic driven model (ModelC).
Table 7.1: Estimated time required for model generation and simulation.
Model A B C D
Segmenting (days) 5 5 4 4
Meshing (days) 1.5 1.5 1 1
Materials (days) 2 2 1 1
Motion data (days) 1 1 1 1
Musculoskeletal modeling (days) 5 5 - 5
Boundary conditions (days) 1 1 0.5 0.5
Converged solution (days) 90 75 4 4
Run time (days) 0.5 0.5 0.45 0.45
Total (days) 106 91 12 17
7.2 FE MODEL RESULTS AND T2/T1ρ RELAXATION TIMES