Towards Patient-Specific Computational Modelling of Articular Cartilage on the Basis of Advanced Multiparametric MRI Techniques

Cartilage degeneration is associated with tissue softening and represents the hallmark change of osteoarthritis. Advanced quantitative Magnetic Resonance Imaging (qMRI) techniques allow the assessment of subtle tissue changes not only of structure and morphology but also of composition. Yet, the relation between qMRI parameters on the one hand and microstructure, composition and the resulting functional tissue properties on the other hand remain to be defined. To this end, a Finite-Element framework was developed based on an anisotropic constitutive model of cartilage informed by sample-specific multiparametric qMRI maps, obtained for eight osteochondral samples on a clinical 3.0 T MRI scanner. For reference, the same samples were subjected to confined compression tests to evaluate stiffness and compressibility. Moreover, the Mankin score as an indicator of histological tissue degeneration was determined. The constitutive model was optimized against the resulting stress responses and informed solely by the sample-specific qMRI parameter maps. Thereby, the biomechanical properties of individual samples could be captured with good-to-excellent accuracy (mean R2 [square of Pearson’s correlation coefficient]: 0.966, range [min, max]: 0.904, 0.993; mean Ω [relative approximated error]: 33%, range [min, max]: 20%, 47%). Thus, advanced qMRI techniques may be complemented by the developed computational model of cartilage to comprehensively evaluate the functional dimension of non-invasively obtained imaging biomarkers. Thereby, cartilage degeneration can be perspectively evaluated in the context of imaging and biomechanics.

In the context of personalized medicine, biomechanical computational modelling becomes ever more relevant in patient-specific care 1 . Modelling-based predictions of biomechanical tissue properties hold great potential for the non-invasive and non-destructive characterization of the tissue status in health and disease; however, these go along with high requirements for simulation and measurement techniques. In particular, the constitutive model has to address significant variations in mechanical properties associated with age, gender, lifestyle as well as disease 2 . In the field of cartilage modeling, patient-specific predictions could significantly facilitate the detection of cartilage degeneration, which is the hallmark change of osteoarthritis (OA). Detection of early changes of particular relevance as the degenerative cascade may at least be slowed (if not even halted) in its progression if early preventive action such as modification of activity level, weight loss, pharmacological chondroprotection or axis-modifying surgery is taken in a timely manner 3 . However, up to now, it is not possible to detect the early, potentially reversible stages of cartilage degeneration using clinical standard imaging modalities 4,5 . Against this background, the non-invasive imaging of cartilage by advanced MRI techniques has made considerable progress over the last decades. Functional MRI techniques (synonymous with quantitative MRI [qMRI]) such as T2, T2* and T1ρ mapping have been developed and validated in a variety of scientific and clinical contexts to characterize extracellular matrix properties of cartilage [6][7][8] and provide measures related to tissue composition and structure 4,9 . T2 relaxation describes the decay of transverse magnetization and is reflective of the tissue's collagen structure and water content. In addition to the tissue-characteristic T2 relaxation, T2* relaxation is governed by additional T2* decay secondary to static magnetic field non-uniformity. In contrast, T1ρ relaxation is determined by measuring the decay of locked transverse magnetization and commonly considered indicative of low-frequency interactions between the tissue's macromolecules and extracellular water. Meanwhile, the exact sensitivity and specificity profiles of both T2* and T1ρ need additional clarification 10,11 . To date, a solid body of evidence has been accumulated indicating the diagnostic potential of such qMRI techniques (excellently reviewed in 6,8,12 ).
Even though previous studies have reported correlations between qMRI parameters on the one hand and structural as well as compositional features of the tissue on the other hand 13,14 , there is clearly no consensus on any specific relations between these parameters. Previously, our group studied correlations between measured qMRI parameter maps and modelled volume fractions to better refine each qMRI parameter's sensitivity and specificity profile 15 . Additionally, quantitative T2 maps have been referenced to loading-induced changes in cartilage composition and structure based on a constitutive cartilage model 16 . However, to the best of our knowledge, multiparametric qMRI maps have not been modelled as a function of the functional properties of the tissue in general and its biomechanical measures in particular. Thus, this study aimed to establish a framework to integrate sample-specific multiparametric qMRI maps into the proposed constitutive material model while subsequently optimizing the model in terms of weighted structural and compositional tissue features as derived from the multiparametric qMRI maps. Therefore, the working hypotheses of the study were that (1) the complex relation between the stress response of the tissue and its qMRI appearance in terms of respective T1, T1ρ, T2 and T2* maps may be translated into a refined constitutive model of cartilage tissue and (2) the functional properties of the tissue can be described by this model following its comprehensive optimization.

Results
Upon histological assessment, all samples were found to be grossly intact (Mankin Grade 0, mean sum score 3.2 ± 0.8).

MRI measurements.
Spatially resolved quantitative T1, T2, T1ρ and T2* maps were obtained for all osteochondral samples. Qualitatively, samples were relatively homogeneous, while uniform changes in signal intensities were found as a function of tissue depth (Fig. 1). Whenever present, focal signal alterations were only slight, and, in any case, adjacent cartilage areas were not affected.
Quantitatively, qMRI parameter values were characterized by considerable standard deviations even though the mean parameter values were comparable. Table 1 presents a detailed overview of the qMRI parameter values. Fig. 2 as a function of time for two representative samples. Considerable variability in the stress responses of the samples can be seen. When keeping the sample strained, larger stress relaxation was observed with larger initial stresses.

Confined compression tests. The nominal stress response is illustrated in
In Fig. 3a nominal stresses are plotted versus the intra-tissue volume changes. Similar to the stress-time response, the stress-induced volume changes were highly variable. Correspondingly, large standard deviations were found when calculating the mean stiffness at ε = 15%: E = 1.53 ± 0.68 MPa (mean ± SD). Sample-specific biomechanical testing details are displayed in Table 1.

Model evaluation.
Relaxed stress calculated by the model is plotted in Fig. 3b as stress versus intra-tissue volume changes. Predictions of the computational model (Fig. 3b) demonstrated grossly similar characteristics (in terms of trend and overall curve characteristics) as compared to the experimental measurements (Fig. 3a).   Table 1. Quantitative characterization of human articular cartilage samples (n = 8) by qMRI parameters and the stiffness at a tissue strain of ε = 15%. The entire sample cross-sectional area was the region-of-interest. Data are given as mean ± standard deviation, while goodness-of-fit measures R 2 and Ω detailing the correspondence between experimentally measured and theoretically modelled data are shown in the last two columns (from right).   Table 1.
The global set of material parameters obtained by the inverse Finite-Element optimization is given in Table 2. By assuming linear relations between the qMRI parameter-derived specific volume fractions φ ξ  x T (x ∈ {1, 1ρ, 2, 2*}) and scalar coefficients ξ w Tx , sample-specific volume fractions Φ ξ for the fluid and collagen contents (ξ ∈ {f, co}) were obtained as a function of normalized sample depth z (see Fig. 4). Towards the sample surface, fluid content was high, while collagen and proteoglycan contents were low. Towards the sample bottom (i.e. the cartilage-bone transition), the differences in volume fractions were less pronounced even though fluid was still the dominant tissue component. Following optimization of the spatial volume fractions φ ξ  z ( ) x T ( Fig. 5a-d), the φ ξ (Tx) values of the qMRI parameters showed close correspondence to the idealized model of the depth-related volume fractions as proposed by Wilson et al. 13 (Fig. 5e). Accordingly, ξ w Tx represent weighting parameters which define each qMRI parameter's contribution to the collagen and fluid volume fractions (please see MRI-based model input section below for more details). Additionally, six constitutive model parameters were obtained by the optimization procedure. Here, k 1 , k 3 and k 3 describe the material non-linearity of the collagenous part, while a 0 and a 1 are associated with the stiffness and compressibility of the non-collagenous matrix. w ∈ [0, 1] denotes a weighting parameter of the alignment of the collagen fibrils towards its preferred direction with the lower limit representing an ideal fibre alignment, while the upper limit expresses an isotropic fiber distribution.

Weighting parameters
Material model parameters  Table 2. Details of the global material parameter set obtained by the inverse Finite-Element optimization of the spatial volume fractions in the framework of the computational model of cartilage. Weighting parameters ( ξ w Tx ) detail each qMRI parameter's contribution to the collagen and fluid volume fractions and are to be read as follows: The fluid content of the tissue is primarily represented by T1 (41%) and T1ρ (42%), while the collagen content is primarily represented by T1 (31%) and T2* (46%). Six material model parameters complement the proposed constitutive model: k 1 , k 3 and k 3 are related to the biomechanical behaviour of collagen fibrils under tension, while a 0 and a 1 are associated with the stiffness and compressibility of the non-collagenous matrix. w describes the concentration of the collagen fibrils in their preferred direction. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The most important finding of this study is that sample-specific and spatially resolved qMRI data may be integrated into a computational model of cartilage to (1) emulate structural and compositional tissue parameters and to (2) reliably capture cartilage functional properties. The computational model of cartilage has been further refined and optimized to integrate imaging (i.e. the qMRI parameter maps) and biomechanical (i.e. its stress response to confined loading) information. www.nature.com/scientificreports www.nature.com/scientificreports/ To the best of our knowledge, this is the first study to bring together imaging and functional tissue parameters within the framework of a computational model. QMRI data were used as input variables in a sample-specific manner, while the remaining parameters were kept constant in efforts to keep the model complexity manageable. In practical terms, the qMRI parameter maps were used to derive measures of structural and compositional tissue properties as a function of sample depth. To this end, the qMRI parameters have been weighted in their respective contributions to the tissue properties on the basis of an idealized cartilage model as proposed by Wilson et al. 17 . Since specific relations between the exact properties of the tissue and the qMRI parameters remain disputed and considerable overlap in specificity has been reported 10,15,16,18-20 the first step was to identify weighting factors for the optimized qMRI-based tissue assessment.
In spite of the thorough and user-independent optimization against functional properties of the tissue, this framework can only provide the starting point to further study the correlations between imaging, compositional and biomechanical parameters. In particular, refined biochemical techniques providing spatially resolved measures of the tissue features such as microspectroscopy or polarized light microscopy (e.g. 21,22 ) should be included in future studies to integrate sample-specific data on exact tissue properties rather than an idealized tissue model that gives tissue constituents as a function of depth. Nonetheless, the Wilson model of cartilage is validated and commonly used to convey depth-related information on tissue properties of cartilage 17,23-25 . Correspondingly, higher contents of extracellular fluid and lower contents of the solid constituents (i.e. collagen and proteoglycans) were found at superficial sample zones, while the opposite was observed for deeper cartilage zones. These compositional features are well in line with published reference data, e.g. by Wilson et al. 23 .
The computational model of cartilage was efficient enough to successfully describe biomechanical properties (e.g. stiffness) as obtained in the subsequent confined compression tests. Herein, both the non-linearity of the stress-relaxation tests as well as the biomechanical quantities determined at equilibrium (both measured and modelled) are reflective of earlier literature findings 26 . Even though the sample size of the present study is limited and the computational model of cartilage is not yet optimized in terms of refined compositional input variables (as outlined above), our results substantiate the fact that advanced qMRI techniques are adequate to determine relevant tissue features (in structure and composition) that determine the biomechanical properties and stress resilience of the tissue and may be used to non-invasively study tissue functionality. Once further optimized, the computational model might be used to complement clinical-standard MRI examinations to provide a detailed representation of the functional properties of the tissue and to identify tissue regions at risk of incipient degenerative changes. This is of particular relevance as articular cartilage is exquisitely sensitive to the mechanical environment and cartilage degeneration in OA is considered as the key pathophysiological result of abnormal mechanics 27,28 . In the clinical context, such patient-specific modelling approaches aim to obtain image-based surrogate parameters of tissue functionality to identify tissue areas at risk Nonetheless, these approaches need further refinement, in particular when translating the findings to the in-vivo and entire-joint configuration.
Although predictions of the proposed model were confirmed to a large extent by the experimental data, there was some discrepancy between modelled and measured datasets. Possible reasons involve a variety of as yet ill-controlled variables: Particular care was taken to standardize storage conditions, yet systematic error may have been introduced by the prolonged a storage in sample non-physiological environment. This may have led to alterations in the extracellular fluid content and biomechanical properties. Here, comparative longitudinal imaging-based assessment of cartilage functionality as a function of loading may be integrated into the model to further refine the input of data. Additionally, the distinct collagen network properties (in terms of orientation and integrity) have not been considered as additional input variables. In view of the ongoing scientific controversy on the imaging correlates of the collagen network 15,29 , future studies should take into account its complex features beyond its mere content, especially when the mechanical behaviour of the entire tissue is of interest. Here, the combination of mechanical stimuli and advanced qMRI techniques may help along since mechanical loading may be applied during imaging to study loading-induced intra-tissue changes to provide additional functional information. This can improve the computational model and its descriptive capacities 11,13,14,30 . Histology was used to confirm gross structural integrity of the cartilage tissue. However, samples may have been affected (in a functionally relevant context) beyond the histologically assessable scope, because OA is a disease that affects the entire joint by triggering catabolic and inflammatory processes in all compartments. As samples were harvested from total knee replacements only, our results remain to be confirmed in truly healthy cartilage tissue, e.g. from organ donor networks or tumor endoprothesis. Moreover, additional research activities have to be aimed at the inclusion of larger sample numbers of variable degeneration (as controlled by histology) to corroborate the potential of the model in predicting functional properties of cartilage in health and disease. Further limitations involve sample size and study setup. In this study, all measured data were used for thorough model parameter optimization. Larger sample sizes, including new samples, need to be included in future studies to assess the model's predictive capabilities. To this end, the sample-specific qMRI data will be the only input data to the model, while the measured respective mechanical responses will be used as benchmark features. By applying k-fold cross validation schemes as in machine learning techniques (e.g. 31 ) the model's predictive capabilities can be assessed. Additionally, the confined compression tests used for biomechanical reference evaluation is unlike the actual weight bearing in vivo, thereby limiting our study's transferability to the in-vivo setting. Against this background, more physiological forms of loading, quite possibly under simultaneous imaging 14,32 , should be applied in future studies.
In conclusion, this study introduces a computational model of cartilage that integrates qMRIN parameter maps as spatially resolved measures of the structural and compositional properties of the tissue to describe its biomechanical properties. On the basis of this model, advanced qMRI techniques may be complemented to comprehensively evaluate the functional dimension of non-invasively obtained imaging biomarkers. Thereby, cartilage (2019) 9:7172 | https://doi.org/10.1038/s41598-019-43389-y www.nature.com/scientificreports www.nature.com/scientificreports/ degeneration (as the hallmark change of OA) may be appreciated in the context of abnormal mechanics and used as a potential target in diagnosing early (and potentially reversible) OA.

Materials and Methods
study design. This study was designed as a prospective, comparative, intra-individual ex-vivo study that aimed to integrate the functional biomechanical properties of cartilage tissue and their qMRI correlates as input variables into the framework of a computational model of the tissue. Prior to this study, local Institutional Review Board approval from the Ethical Committee of RWTH Aachen University, Germany (AZ-EK157/13) was obtained. Only after individual oral and written informed patient consent was the material that had been collected intraoperatively included in the present study. Moreover, all consecutive experiments were performed in accordance with relevant guidelines and regulations.
Cartilage sample preparation. Human articular osteochondral samples were prepared as in earlier studies 11,14,20 . Briefly, macroscopically intact osteochondral samples were obtained from eight consecutive patients undergoing total knee replacement at our institution from 11/2017 to 03/2018 (3 male, 5 female; age 64.5 ± 12.3 years [range, 40-77 years]). By obtaining one sample from one patient sample pooling was avoided. The osteochondral samples included in this study were harvested from patients with primary osteoarthritis. Therefore, all forms of secondary OA or other bone and joint disorders were excluded. After its sterile excision during surgery, the material was collected in Dulbecco's modified Eagle's medium (DMEM, Gibco-BRL, Gaithersburg, MD, USA) with a set of standard antibiotics added (i.e. 100 U/ml penicillin [Gibco], 100 μg/ml gentamycin [Gibco] and 1.25 U/ml amphotericin B [Gibco]). The osteochondral samples were then prepared according to the standard procedure as follows: First, for the sake of topoanatomic consistency, samples were harvested from the lateral femoral condyles only. Second, samples were evaluated macroscopically according to the Outerbridge classification 33 and only grossly intact samples were included (i.e. Outerbridge grades 0 and 1). Structural integrity was subsequently confirmed by means of histology (see below). Third, samples were cut to standard square shape (20 × 20 mm [width × length]) and the subchondral lamella was preserved, while all cancellous bone was removed. Samples had a mean thickness of 4.201 mm ±1.01 mm (Mean ± Standard Deviation). Using a rongeur, two notches were created at opposing sample sides to define the mid-sagittal plane for future reference.

MRI measurements, data acquisition and analysis. After sample preparation, MR imaging exam-
inations were performed on a per-sample basis using a clinical standard 3.0 T MRI scanner (Achieva, Philips, Best, The Netherlands). As before 11,30 , samples were placed in a transparent container fully immersed in DMEM solution with additives. Samples were imaged using a modified single-channel prostate coil (BPX-30 disposable endorectal coil, Medrad/ Bayer, Leverkusen, Germany) 14,20,34 . Particular attention was paid to position the samples at the iso-center of the coil while aligning the mid-sagittal plane along and the surface parallel to the main magnetic field B 0 . Prior to scanning, B 0 inhomogeneities were excluded using B 0 mapping. After scout views, proton-density weighted sequences were acquired in the axial, coronal and sagittal planes oriented perpendicular to each other (Table 3). On the basis of the axial views, the sagittal imaging section was guided along the mid-sagittal plane to generate a centrally bisecting plane through the sample. Afterwards, T1, T1ρ, T2, and T2* sequences were acquired with the sequence parameters detailed in Table 3. MR imaging was performed at room temperature monitored before and after the measurements (19.5 ± 0.7 °C). Once the data acquisition was completed, the MR raw data including time constants for each pixel were loaded into Matlab R2016a software (Natick, MA, USA). Then, spatially resolved parameter maps were generated by means of predefined and customized fitting routines on a per-pixel basis. Individual pixel values were determined with T1, T1ρ, T2 and T2* relaxation times calculated as follows T1, T1ρ, T2 and T2* were the target relaxation times to be quantified on a per-pixel basis. Here, T E is the echo time, T SL the duration of the spin lock pulses, T R the repetition time and T I the inversion recovery time (i.e. the time delay between the initial inversion recovery pulse and the read out), while A and B are the signal pre-factor and offset, respectively, accounting for proton density and background noise. R 2 statistics adjusted to the number of degrees of freedom was used to check the quality of the fits. For T2 and T2*, only pixel values of expected echo times (T E ≤ 60 ms) were included to reduce the potential of mis-fitting. Sample segmentation was performed manually on the proton density-weighted morphologic images by choosing pixels that safely lay within the tissue.
www.nature.com/scientificreports www.nature.com/scientificreports/ Boundary pixels were excluded to avoid partial volume effects. Segmentation of sample outlines was subsequently validated against the parameter overlays.

Confined compression tests.
Within 24 h after the MR measurements, confined compression tests were conducted using a custom-made compression device. In line with literature data 26,35 , the osteochondral samples were placed in an impermeable confining chamber matching the dimensions of the samples (with a diameter of d = 8 mm). Porous stainless steel filters were placed above and underneath the samples to allow for fluid outflow. The samples were compressed by a hollow piston which axially moved down the upper filter (Fig. 6). The entire set-up was placed in a container filled with phosphate-buffered saline (PBS, Gibco) and mounted on a universal testing machine (Zwick Roell Z010, Zwick GmbH, Ulm, Germany).
Prior to the tests, a tare load of 1 N was applied to the osteochondral sample by lowering the piston at a rate of 0.16 mm/min to maintain standardized interaction of the interfaces throughout the measurements 26 . The resulting piston position was used to determine sample thickness, while equilibration was achieved by holding the piston in place for 15 min. Then, the osteochondral samples were subjected to a sequence of 20 ramped compressions with a strain-rate of 1%/min to a total strain of 20%. After each applied loading step, relaxation phases variable in duration (as depicted in Fig. 2) were applied to allow for sufficient equilibration 26 .
Histological analyses. After biomechanical tests, samples underwent histological processing by simultaneous decalcification and fixation (Ossa fixona, Diagonal, Münster, Germany), sectioning along the mid-sagittal plane and embedding in paraffin. From the mid-sagittal plane, 5-μm-thick sections were cut, stained with hematoxylin/eosin and Safranin O according to standard protocols and visualized using a Leica light microscope (model DM LM/P, Wetzlar, Germany). Histological grading was performed in a blinded manner using the Mankin classification 36 . Based on the quantitative assessment of tissue structure, cellularity, proteoglycan staining intensity and tidemark integrity, the Mankin sum score is a representative measure of tissue degeneration. Ranging from 0 to 14, a Mankin sum score of 0 is indicative of no histological signs of degeneration. Mankin sum scores may be grouped into distinct Mankin grades; i.e. Mankin grade 0 indicates structurally grossly intact cartilage (Mankin sum scores 0-4) 37 .    www.nature.com/scientificreports www.nature.com/scientificreports/ MRI-based model input. The qMRI parameter maps were used as exclusive sample-specific input values within the finite element (FE) code. Material constants as well as the weighted qMRI parameters (see Table 2) were obtained in an inverse FE manner with an optimization algorithm updating the parameters for all samples simultaneously and thereby enforcing a global set of parameters.

Sequence Type Parameters
Cartilage is classically considered as a bi-phasic material with a solid phase hydrated by interstitial fluid. The fluid-filled solid extracellular matrix primarily consists of the collagen (CO) fibril and proteoglycan (PG) fractions, in the following represented with respect to volume by φ co and φ pg , respectively 17 . The interstitial fluid and solid volume fractions are denoted by φ f and φ s , respectively. They are related by φ s = 1 − φ f . The volume fractions are expressed as a function of the qMRI parameters Tx, which satisfy the following condition: Functional dependencies between qMRI parameters and idealized cartilage volume fractions as characterized in previous studies 15,16 provided the basis for the qMRI-informed computational model of cartilage. More specifically, relations between the qMRI parameters and the fractional composition of an idealized cartilage model (as defined by Wilson et al. 23 ) were created in a pixel-wise and depth-related manner for the fluid and CO fractions (φ f and φ co , respectively). Note that this model of cartilage composition defines the mean content of each cartilage constituent as a function of tissue depth. These datasets served as the basis for the non-linear regression analysis with an exponential dependency where a, b and c are fitting coefficients. Based on earlier data 15,16 , these coefficients were determined by optimizing the contribution of each cartilage constituent to the qMRI parameter maps. Accordingly, the inverted generalized form is denoted by φ ξ (Tx), where the function is continuously extended by the constant parameters l Tx,ξ and u Tx,ξ at the lower Thus, one arrives at four distinct qMRI parameter specific volume fractions for CO and fluid, which were fitted against the pixel-wise resolved qMRI parameters. To this end, a one-dimensional phenomenological representation over the normalized depth coordinate z was employed as As all cartilage tissue phases are considered incompressible, any volume decrease is assumed to be solely due to fluid outflow 39,40 . Hence, after complete fluid loss the fluid volume fraction is approximately 0 and the tissue volume approaches its solid constituent fraction (also referred to as compaction point at which J → Φ s ) 41,42 . For a bi-phasic material the Helmholtz free energy (cf. 43 ) can be given as where the first term on the right hand side (Ψ C ( )) reflects the isochoric free energy resulting from the deformation of the solid parts, while the second term (Ψ π (J, μ)) describes the chemo-mechanical interactions. The third term (μC) is determined by the chemical potential of the solvent μ and the molar solvent concentration = Φ C J V f m with the molar volume V m (cf. 43,44 ). Hence, C governs the fluid flux and is consequently neglected in the following, since the modelling is restricted to the relaxed (i.e. time-independent) configuration. Accordingly, the Cauchy stress tensor is given as σ = ∂ Ψ ∂ + ∂Ψ ∂ 2/3 1 1 while S represents the isochoric contribution of the second Piola-Kirchhoff stress tensor. Furthermore, π denotes the osmotic multiplier and can be given by an empiric expression (cf. 43 ) as (11) 2 is a phenomenological relation enforcing a stronger correlation with the sample-specific qMRI data. a 0 and a 1 are material parameters, while π 0 ensures a stress-free reference configuration. As the extracellular matrix is composed of the solid cartilage phase (i.e. CO and PG) 23,39,40 , its free energy can be given as with the parameter a 0 which enforces a stress-free reference configuration, see (11) 1 . = I C tr 1 denotes the first isochoric strain invariant. The fibril network is modelled by a polyconvex expression which can be given as (cf. 46 captures the J-shape response to tension typical for CO fibres, while the latter term −   ( ) describes the fibre contribution to the tissue response to compression due to the tube contraction effect 48,49 . k 1 , k 2 and k 3 denote material constants. The local CO fibril architecture is largely responsible for cartilage anisotropy. Hence, in the constitutive model, cartilage anisotropy is captured by weighted structural tensors L i and associated structural invariants where w i (i = 1, 2, ..., m) denote scalar weighting parameters associated with the preferred fibre families whose directions are specified by unit vectors m i . I denotes the identity tensor. Finally, in view of (10) 3 and (12), the isochoric part of the second Piola-Kirchhoff stress tensor can be given as Numerical implementation and parameter estimation. For the parameter identification, the computational model of cartilage proposed in the previous section was implemented as a user-subroutine UMAT (user material routine) into the commercially available software package Abaqus FEA (v6.17, Simulia Corp., Providence, USA). To emulate the idealized Arcade model of Benninghoff 51 , the preferred fibril directions θ fib (z) were implemented in a rotationally symmetric manner with eight equiangular fibre families as before 16 . The weighted fractional relations (7) were mapped onto the sample-specific geometry. To exclude the possibility of negative volume fractions as a result of the saturation condition (5), the condition of minimally positive PG fraction values was enforced by the tolerance value Δtol = 1 × 10 −6 . Algorithm 1 presents a detailed overview of the code used for FE implementation of the computational model of cartilage. In the confined compression test the sample-specific geometry was implemented on the basis of its mean thickness, while the piston was modeled as a rigid body.
The parameter fitting was performed using a constrained optimization in Matlab while satisfying condition (8). The minimal error between the FE simulations and actual experiments for the complete set of samples (n = 8) served as the objective function. An overview of the optimization framework can be found in Algorithm 2.
The accuracy of the simulations as compared to the experimental measurements was quantified by determining the square of Pearson's correlation coefficient (i.e. R 2 as in 52 ) and the relative approximated error defined by Ω = ⋅ − y f y 100 ( ) / , where y and f denote vectors of experimental and computational values, respectively.