Elastic, Dynamic Viscoelastic and Model-Derived Fibril-Reinforced Poroelastic Mechanical Properties of Normal and Osteoarthritic Human Femoral Condyle Cartilage

Osteoarthritis (OA) degrades articular cartilage and weakens its function. Modern fibril-reinforced poroelastic (FRPE) computational models can distinguish the mechanical properties of main cartilage constituents, namely collagen, proteoglycans, and fluid, thus, they can precisely characterize the complex mechanical behavior of the tissue. However, these properties are not known for human femoral condyle cartilage. Therefore, we aimed to characterize them from human subjects undergoing knee replacement and from deceased donors without known OA. Multi-step stress-relaxation measurements coupled with sample-specific finite element analyses were conducted to obtain the FRPE material properties. Samples were graded using OARSI scoring to determine the severity of histopathological cartilage degradation. The results suggest that alterations in the FRPE properties are not evident in the moderate stages of cartilage degradation (OARSI 2-3) as compared with normal tissue (OARSI 0-1). Drastic deterioration of the FRPE properties was observed in severely degraded cartilage (OARSI 4). We also found that the FRPE properties of femoral condyle cartilage related to the collagen network (initial fibril-network modulus) and proteoglycan matrix (non-fibrillar matrix modulus) were greater compared to tibial and patellar cartilage in OA. These findings may inform cartilage tissue-engineering efforts and help to improve the accuracy of cartilage representations in computational knee joint models. Supplementary Information The online version contains supplementary material available at 10.1007/s10439-021-02838-4.


Details on sample preparation
Thirty-five osteochondral samples from medial and lateral femoral condyle of total knee replacement (TKR) patients due to medial compartment tibiofemoral OA (number of samples = 16, number of subjects = 14, age 64.5 ± 7.9 years, age range 53 79 years, 6 men and 8 women) and healthy donors without known clinical knee OA or rheumatoid arthritis (number of samples = 19, number of subjects = 10, age 53.4 ± 18.8 years, age range 18 77 years, 5 men and 5 women) were obtained by a hole saw (diameter = 4 mm). Initially, our sample set consisted of 50 samples but most of the medial condyle samples from TKR patients and one sample from healthy donors were not measurable in biomechanical indentation geometry due to full-thickness cartilage loss, thus they were not included in the study.
Surgeries were conducted in Trelleborg Hospital, Trelleborg, Sweden while donor samples were obtained from Skåne University Hospital, Lund, Sweden. All samples were obtained within 48 h postmortem and frozen at -80 o C within 2 hours of extraction in phosphate buffered saline (PBS). Afterwards, the samples were transferred from Sweden to Finland, while kept frozen in dry ice (~-78 0 C). They were then kept at -23 o C until conducting experiments. The storage of samples at this temperature should not cause changes in composition 25,28 or mechanical properties of the tissue 27 .

Details on mechanical indentation testing
The osteochondral sample in PBS was thawed at room temperature for 15 minutes before starting the measurements. Based on literature 27 , this was considered enough to thaw the sample. Cartilage thickness from each osteochondral sample was measured from the surface to the cartilage-bone interface using an optical stereomicroscope (9 1.6, Zeiss, STEMI, SV8, Germany). The measurements were conducted from 4 different quarters on the sample side. The mean value was then calculated to represent the sample thickness. We then flattened the bone surface using sandpaper (Mirox P80,Mirka Oy,Uusikaarlepyy,Finland) to reduce roughness and increase adherence of bone to the measurement chamber. The bone end was then glued (Loctite Precision, Henkel AG, Düsseldorf, Germany) on the metallic bottom of a custom-made transparent acrylic chamber, which was filled with PBS.
A custom-made high-precision linear servo-motorized indentation device (Newport PM500-C Precision Motion Controller, Newport PM1A1798 Actuator, Irvine, CA, USA) equipped with a 1000 g load cell (Honeywell afterwaModel 31/AL311BL, Columbus, OH, USA) and a custom-made cylindrical plane-ended metal indenter (diameter = 0.73 mm) were used to measure the mechanical response of the samples. This indenter diameter is small enough to minimize effects from the sample edges 26 .
Prior to starting the measurement, a pre-stress of 12.5 kPa was applied to establish a proper sample-indenter contact for consistent and repeatable measurements, similarly as in previous studies 13,19 . We allowed the sample to equilibrate for about 15 minutes 3,9 and checked that pre-stress remained at 12.5 kPa. We then set the zero-strain level (reference strain). Subsequently, a 4-step stress-relaxation protocol was applied. Each step included the application of a 5% compressive strain (of the remaining cartilage thickness) followed by 15 minutes of relaxation 19 . In preliminary tests, we observed that this relaxation time is enough to reach equilibrium. We applied strain rate of 100%/s to ensure maximum fluid pressurization (and collagen fibers involvement) in the tissue. The same strain rate has also been used in previous studies 19,21 . A recent review article suggested that testing at higher strain rates is dominated by initial fluid pressurization, whereas slower strain rates (e.g., 0.1%/s) allow for some energy dissipation by fluid flow 24 . Moreover, strain-rates of this magnitude are known to occur during daily activities 2 . The multi-step protocol was chosen to investigate the nonlinear properties from the peak forces (i.e. nonlinear fibril network and instantaneous modulus) and relaxation phases (i.e. nonlinear strain-dependent permeability) as well as the linear stress-strain response from the equilibrium points. Some in-vivo studies have reported nominal cartilage strain of up to 20% during gait cycle 1,18,22 . Moreover, previous studies have suggested that ~20% total strain does not cause plastic deformation or damage to cartilage 13,19 . After the stress-relaxation measurement (the sample at the final strain of the stress-relaxation experiment), a dynamic sinusoidal test was conducted using a 2% strain amplitude (of the remaining thickness) with frequencies of 0.005, 0.05, 0.1, 0.25, 0.5, 0.625, 0.833 and 1 Hz (four cycles for each frequency). Based on earlier studies 16,29 , this frequency regime was assumed to demonstrate frequency-dependent changes in the dynamic modulus and phase difference. The 2% strain level provided us with a measurable response in the cartilage tissue, while keeping the total strain level within the safe region 13,19 .

Details on elastic and dynamic viscoelastic properties
The equilibrium modulus of cartilage was calculated from a slope of a linear leastsquares fit to the stress-strain points at the equilibrium (Figure 1, middle column). The modulus was corrected using the Hayes correction factor 6 as follows: where -dimensional constant .e. a/h)). We set 11 . To obtain the dynamic moduli at different frequencies, the stress and strain amplitudes were calculated from each cycle and averaged over the 4 measured consecutive cycles. We also corrected the dynamic moduli using the Hayes correction factor 6 . Here, we assumed the sample to be incompressible during the dynamic loading (as fluid is pressurized due to high loading rate), thus 10 . The Fourier transform was used to extract the frequency content of dynamic data. The phase angle was extracted at a frequency at which the power amplitude peaks. The phase difference was calculated by subtracting the displacement and force phase angles. The instantaneous modulus was calculated from the peak stress values at each step ( Figure 1, middle column). Thus, we obtained data points for the instantaneous modulus as a function of strain.
Accordingly, a linear least-squares line was fitted to the data points of instantaneous modulus vs. applied strain and used to determine the initial instantaneous modulus (intercept) and strain-dependent instantaneous modulus (slope). These moduli were also corrected using the Hayes correction factor 6 with 10 as the sample was assumed to be incompressible due to the high loading rate.

Material model
Cartilage tissue was modeled using the FRPE material model, in which articular cartilage is composed of an elastic fibrillar matrix (representing the collagen fiber network) and a porous hyperelastic non-fibrillar matrix (representing the PG matrix), filled with fluid.
The collagen network was modeled with 4 organized collagen fibrils (primary fibrils) and 13 randomly oriented fibrils (secondary fibrils) 19,32 . The density ratio of the primary fibrils to the secondary fibrils was fixed to 12.16 8,19 . The compressive stress-strain behavior of the fibrils was set to zero, while in tension the following non-linear stressstrain behavior was assumed: where and are stress and strain of the fibril, is the initial fibril network modulus and is the strain-dependent fibril network modulus 17 . Furthermore, the non-fibrillar matrix was modeled using the Neo-Hookean hyperelastic material model. The below formulation for the non-fibrillar matrix stress was used: where is the stress tensor of the non-fibrillar matrix, and are the shear and bulk moduli of the non-fibrillar matrix, respectively, F is the deformation gradient tensor, J is the determinant of the F and I is the unit tensor 31 . The bulk and shear moduli of the non-( ) and Poisson ratio ( , based on 12,19 ) of the non-fibrillar matrix: Regarding the fluid flow behavior, 7 was recruited to describe the fluid flow inside the porous matrix as follows: where is the fluid flow flux, is the (hydraulic) permeability of the material and is the (fluid) pressure gradient.
which is true in most biological tissues 5 . The deformation in the porous material causes a change in the void ratio (the proportion of the fluid volume to the solid volume), and consequently, changes in the permeability 30 , which is formulated in the following way: where and are the current and initial values for the permeability, and and are the current and initial values for the void ratio, respectively. M is a constant describing the void-ratio (or deformation) -dependency of permeability 5,30,32 . The initial void ratio which is the ratio of the fluid volume to the solid volume was set to 3, based on literature 14 .

Model construction and boundary conditions
Sample-specific axisymmetric models were built in Abaqus (V6.14, Dassault Systèmes Simulia Corp., Providence, RI). The samples were meshed using linear axisymmetric pore pressure continuum elements (element type CAX4P). A mesh convergence test was conducted to ensure a proper mesh size. The structure and composition were assumed homogenous (i.e. collagen fiber network orientation parallel to the cartilage surface, homogeneous PG and collagen contents, and constant void ratio) to obtain material properties independent from the composition and structure of the tissue, similarly as done before 4,12,20 . The compression of the indenter was simulated by a displacement boundary condition on the cartilage surface for computational efficacy.
The contact between the lateral edge of the indenter (an analytical rigid surface) and cartilage surface was modeled using a frictionless hard contact (both in normal and tangential direction; during contact the separation of surfaces was allowed only in tangential directions) to prevent folding of the cartilage mesh. Free draining was allowed from non-contacting surfaces (i.e. pore pressure = 0), while cartilage-indenter contact was modeled impermeable. The axisymmetric boundary condition was applied on the symmetry axis of the sample, i.e., lateral displacement of the nodes at the symmetry axis was fixed and fluid was not allowed to flow through the symmetry axis 19 .
Axial and lateral displacements of the bottom nodes of cartilage were fixed because the subchondral bone was considered rigid. The fluid flow was restricted through the cartilage-bone interface. The models were solved using the soils consolidation analysis in Abaqus.
The following material parameters of the FRPE model ( ) were obtained by fitting the force-time response of the second and third steps of the model to the corresponding steps of the experimental stress-relaxation test. The first stress-relaxation step was simulated in the model but not included in the fitting procedure as this step was considered not to be always in a proper contact 19 . In addition, if more than two steps are used in an optimization procedure, the high level of nonlinearities resulting from inherent inhomogeneities of cartilage may become dominant, thus leading to a poor fit quality with the current material model. The data fitting was conducted using the Nelder-Mead simplex algorithm (

Details on statistical analyses
We analyzed the dependent variables (i.e. the FRPE material parameters ( ) as well as elastic and dynamic viscoelastic material properties) to compare mean values between the groups. A linear mixed-effects (LME) model was used. This statistical model accounts for the dependency of the samples obtained from the same subject (donor or patient). Subjects in each group were set as a mixed (random) effect while the OA progression group (i.e. normal or moderate OA), the compartment (i.e. medial or lateral) and the interaction of the variables were set as fixed variables. The LME models were age and BMI adjusted for group-wise comparisons of femoral condyle cartilage samples.
Regarding the dynamic moduli measured as a function of frequency, we analyzed it using linear splines, because of the non-linear association between frequency and the outcomes (dynamic moduli). The LME model with dynamic moduli as outcomes, and the three linear splines (with knots at 0.05 Hz and 0.5), compartment (medial or lateral) and OA status (moderate OA or normal) as independent variables was used. We included the subject and compartment as random effects and allowed for random slopes for the first spline (as there was no relevant variability in the slopes of the other two splines, see supplementary Figure S3). We assumed an unstructured covariance od for the calculation of degrees of freedom.
For the analysis of phase differences, we treated the frequency as a categorical variable with frequencies of 0.25 Hz and higher as one group, based on inspection of individual trajectories (supplementary Figure S5). We fitted a linear mixed model with compartment, group, frequency, and their interactions as independent variables. We included random intercepts for the subject, compartment and random slopes. We allowed for heteroskedastic errors as the variability was higher at a frequency of 0.005 Hz. We assumed an for the calculation of degrees of freedom. The analyses were repeated with additional adjustments for age and BMI.
Moreover, an age-adjusted LME model was used to compare the FRPE properties of femoral condyle samples to those of other human cartilages (i.e. femoral condyle vs.

Detailed analyses of fibril-reinforced poroelastic (FRPE) parameters
Descriptive data on FRPE material parameters for each group is shown in Figure S1 (each color represents samples from the same knee, except red color, which represents subjects contributing with only one sample). Moreover, descriptive data on FRPE material parameters for each OARSI grade are shown in Table S1. Figure S1: Descriptive data of FRPE parameters grouped according to health status in different locations (medial vs. lateral). Each color represents samples from the same knee, except red color, which represents subjects contributing with only one sample. The results of linear mixed model analyses comparing FRPE parameters between normal and moderate OA groups at medial and lateral femoral condyle cartilage are shown in Table S2.
The differences (with 95%CIs) in fibril-reinforced poro(visco-)elastic properties of femoral condyle cartilage and other sites (tibial, patellar and hip joint cartilage) in different OA stages can be seen in Table S3.

Detailed analyses of equilibrium, initial and strain-dependent instantaneous moduli
Descriptive data on equilibrium, initial and strain-dependent instantaneous moduli for each group shown in figure S2 (each color represents samples from the same knee, except red color, which represents those subject contributing with only one sample).
Moreover, descriptive data on FRPE material parameters for each OARSI grade are shown in Table S4. Figure S2: Equilibrium, initial and strain-dependent instantaneous moduli grouped according to health status (normal and moderate OA) in different locations (medial vs. lateral) The results of linear mixed model analyses comparing elastic parameters between normal and moderate OA groups at medial and lateral femoral condyle cartilage are shown in Table S5.

Detailed analyses of dynamic moduli and phase differences as a function of frequency
The raw dynamic moduli as a function of frequency for two OA status (i.e. normal and moderate OA) groups and locations (i.e. medial and lateral) are shown in figure S3.
Moreover, descriptive data for Elastic dynamic moduli measured at different frequencies are shown for each OARSI grade in Table S6.  The dynamic moduli were analyzed using three linear splines (with knots at 0.05 Hz and 0.5 Hz). The results of the linear mixed model for dynamic modulus are presented in Table S7. Moderate OA vs normal, medial compartment -1.8 (-4.5, 0.9) -1.7 (-4.6, 1.2) Moderate OA vs normal, lateral compartment -1.1 (-3.0, 0.8) -0.8 (-2.9, 1. 3) The statistical model considers the person-specific effects (i.e. which samples are obtained from the same knees). To further clarify it, the dynamic moduli for each subject as a function of frequency is shown in figure S5. As can be seen, samples from  The raw phase differences as a function of frequency for two groups (normal and moderate OA) groups and compartments (medial and lateral) are shown in figure 3. Moreover, descriptive data for dynamic viscoelastic properties (i.e. phase difference) measured at different frequencies are shown for each OARSI grade in Table S8.