Maximum shear strain-based algorithm can predict proteoglycan loss in damaged articular cartilage

Post-traumatic osteoarthritis (PTOA) is a common disease, where the mechanical integrity of articular cartilage is compromised. PTOA can be a result of chondral defects formed due to injurious loading. One of the first changes around defects is proteoglycan depletion. Since there are no methods to restore injured cartilage fully back to its healthy state, preventing the onset and progression of the disease is advisable. However, this is problematic if the disease progression cannot be predicted. Thus, we developed an algorithm to predict proteoglycan loss of injured cartilage by decreasing the fixed charge density (FCD) concentration. We tested several mechanisms based on the local strains or stresses in the tissue for the FCD loss. By choosing the degeneration threshold suggested for inducing chondrocyte apoptosis and cartilage matrix damage, the algorithm driven by the maximum shear strain showed the most substantial FCD losses around the lesion. This is consistent with experimental findings in the literature. We also observed that by using coordinate system-independent strain measures and selecting the degeneration threshold in an ad hoc manner, all the resulting FCD distributions would appear qualitatively similar, i.e., the greatest FCD losses are found at the tissue adjacent to the lesion. The proposed strain-based FCD degeneration algorithm shows a great potential for predicting the progression of PTOA via biomechanical stimuli. This could allow identification of high-risk defects with an increased risk of PTOA progression.


Introduction
Articular cartilage is an avascular soft tissue which reduces the friction between articulating bone surfaces and transmits loads to the underlying subchondral bone. The loadcarrying function of the biphasic cartilage stems from the interplay of the porous solid and mobile liquid phase (Chung and Mansour 2015). From the mechanical point of view, the most important solid constituents in the extracellular matrix (ECM) of cartilage are the proteoglycans (PGs) and the wellorganized arcade-shaped collagen fibrils (Buckwalter and Mankin 1998;Hosseini et al. 2014). The hydrophilic glycosaminoglycans (GAGs) within the PGs create a negative fixed charge density (FCD) in the tissue and attract fluid through osmotic pressure (Urban et al. 1979;Julkunen et al. 2013;Hosseini et al. 2014), thus causing an important swelling behavior of cartilage (Urban et al. 1979;Julkunen et al. 2013). Additionally, they stabilize the ECM and improve tissue resilience and contribute to its durability (Buckwalter and Mankin 1998;Heijink et al. 2012). On the other hand, the collagen network resists tissue swelling and deformation (Basser et al. 1998;Julkunen et al. 2013) and retains the PGs within the tissue ). In the superficial zone of cartilage, the collagen fibrils are oriented parallel to the tissue surface and provide shear resistance for the tissue (Zhu et al. 1993;Motavalli et al. 2014). Together, these components provide the compressive stiffness for cartilage and an ability to carry repetitive external loads in the complex mechanical environment of the knee joint (Buckwalter and 1 3 Mankin 1998; Julkunen et al. 2013). However, an adverse biomechanical environment caused by cartilage lesions and focally elevated contact stresses undermine the function of cartilage. Such detrimental environment formed after a traumatic impact may ultimately lead to a loss of mechanical integrity of the tissue and development of post-traumatic osteoarthritis (PTOA) (Heijink et al. 2012;Mohanraj et al. 2014;Lieberthal et al. 2015). Despite extensive research, various mechanisms behind the development of PTOA are not fully understood Heijink et al. 2012).
Defects penetrating into the cartilage (Hosseini et al. 2014) as well as collagen damage (Stoop et al. 2001;Wilson et al. 2006) cause tissue roughening, thus leading to the loss of its mechanical sustainability (Hosseini et al. 2014). The roughened tissue with increased porosity (Sauerland et al. 2003), higher permeability (Sauerland et al. 2003;Julkunen et al. 2013;Hosseini et al. 2014;Mononen et al. 2018) and damaged collagen network cannot retain PGs similarly to a healthy tissue (McCormack and Mansour 1998;Párraga Quiroga et al. 2017a). Therefore, the high interstitial fluid velocities due to high fluid pressure during compressive loading may lead to PG leakage out of the cartilage through lesion surfaces (Sauerland et al. 2003). Furthermore, the strain levels can increase in the vicinity of cartilage cracks (Gratz et al. 2008;Venäläinen et al. 2016). Such elevated strain levels could alter the metabolic function of cells (Kurz et al. 2001;Gratz et al. 2009) and exceed the limits of chondrocyte (cartilage cell) death (D'Lima et al. 2001;Gratz et al. 2008;Venäläinen et al. 2016;Mononen et al. 2018), collagen cross-link failure (Wilson et al. 2006;Venäläinen et al. 2016) and PG degradation (Hosseini et al. 2014;Venäläinen et al. 2016;Mononen et al. 2018). Local cell apoptosis impairs tissue maintenance (Gratz et al. 2008), for example via reduced PG synthesis rate (Jeffrey et al. 1997). On top of the mechanical factors, the biochemical environment (e.g., inflammation, altered synthesis of cartilage constituents) in the tissue plays a crucial role in cartilage deterioration and PTOA development (Buckwalter and Mankin 1998;Alexopoulos et al. 2005;Heijink et al. 2012;Akkiraju and Nohe 2015;Lieberthal et al. 2015). Altogether, these processes could lead to a local PG depletion and subsequent decrease in the FCD content near the lesion, which is an early sign of osteoarthritis (Heijink et al. 2012;Julkunen et al. 2013;Hosseini et al. 2014).
Since the self-repair ability of cartilage is limited (Akkiraju and Nohe 2015;Makris et al. 2015;Párraga Quiroga et al. 2017a), the effects of tissue defects may accumulate over time and lead to PTOA progression. Biomechanical modeling is an effective way to study the effects of different factors contributing to tissue deterioration as well as to predict failure locations and damage progression. For example, Hosseini et al. (2014) and Párraga Quiroga et al. (2017a) have previously implemented a biomechanical damage progression model based on a strain threshold for degeneration of a PG-rich ground substance. In Stender et al. (2015), the molecular density of the GAGs (i.e., the amount of PGs) was controlled by the volumetric deformation. All in all, strain/deformation has been frequently associated with the PG loss both in experimental (Quinn et al. 2001;Li et al. 2013;Song et al. 2014) and in computational studies (Wilson et al. 2006;Hosseini et al. 2014;Stender et al. 2015;Gardiner et al. 2016;Mononen et al. 2018). However, additional insights into the effects of different strain measures on the degeneration of PGs in PTOA are needed. A robust degeneration algorithm could be a part of a computational framework providing prediction of disease progression and determining high-risk defects [mainly due to the location and size (Peña et al. 2007;Venäläinen et al. 2016)], and ultimately helping in planning of intervention strategies, such as lifestyle change, rehabilitation and surgery.
In order to understand better the effects of mechanical stimuli and the changes that appear around cartilage lesions, the objective of this study is to develop a PG degeneration algorithm which could predict the effect of cartilage defects on the FCD distribution as a function of internal tissue strains/stresses. While the PG degeneration has previously been modeled, for example, by scaling the equilibrium properties of the tissue [e.g., non-fibrillar matrix modulus (Hosseini et al. 2014;Mononen et al. 2018)], relating the FCD concentration directly to the PG degeneration in the proposed algorithm is a novel approach. Based on our previous study (Venäläinen et al. 2016) and degeneration thresholds determined and approximated from the literature (D'Lima et al. 2001;Wilson et al. 2006;Danso et al. 2014;Hosseini et al. 2014;Venäläinen et al. 2016;Párraga Quiroga et al. 2017a;Mononen et al. 2018), we hypothesized that an increase in strain levels near a cartilage injury will lead to a localized decrease in the FCD content in the proximity of the injury. In the algorithm, the FCD loss is controlled by one of the following factors: (1) axial strain, (2) shear strain, (3) octahedral shear strain, (4) maximum shear strain and (5) minimum principal strain. Moreover, we evaluated the effect of (6) maximum principal stress with appropriate degeneration thresholds (Williamson et al. 2003;Eleswarapu et al. 2011;Danso et al. 2014;Mononen et al. 2016;Venäläinen et al. 2016) as it has been suggested to predict subject-specific progression of cartilage degeneration Liukkonen et al. 2017). In the present study, these strain and stress measures are applied and tested for the FCD depletion. Only the FCD loss was considered, since substantial reduction in the PG content has been observed after a trauma (Fick et al. 2016), as opposed to a negligible change in the collagen content (Nelson et al. 2006). Potential collagen fibrillation or disorganization was neglected as the collagen orientation angle was not suggested to alter in the superficial and transitional zones near the cartilage injury (Rolauffs et al. 2010;Tanska et al. 2018). The loss of the FCD via fluid flow or diffusion was not directly modeled, since the literature is lacking the corresponding parameter values. The study hypothesis is investigated, and the performance of the algorithm is evaluated, with finite element (FE) models of intact and damaged 3D cylindrical cartilage explants in unconfined compression.

Material model
Articular cartilage in the 3D model was considered as a fibril-reinforced poroviscoelastic material with Donnan osmotic swelling (FRPVES) (Wilson et al. 2005b). The material model assumes that the cartilage has a biphasic nature with interacting fluid and porous solid matrices. The solid compartment is composed of a viscoelastic fibrillar part, i.e., the collagen network, and a swelling non-fibrillar part comprising primarily of PGs. Depth-dependent mechanical behavior (Szarko and Xia 2012) and tension-compression nonlinearity (Chen et al. 2016) of cartilage are encompassed by the FRPVES model. For the non-fibrillar part, the strain energy density function W nf was defined based on a hyperelastic compressible neo-Hookean solid material  where = T is the right Cauchy-Green deformation tensor and is the deformation gradient tensor. K nf and G nf are the bulk and shear moduli of the non-fibrillar matrix, respectively, and they are defined as where E nf and nf are the Young's modulus and Poisson's ratio of the non-fibrillar matrix, respectively. Thus, the Cauchy stress tensor in the non-fibrillar matrix becomes (Wilson et al. 2005b Fig. 1 Schematic diagram of a viscoelastic collagen fibril. A linear elastic spring (initial fibril network modulus E 0 ) connected in parallel with a nonlinear spring-dashpot system. In the spring-dashpot system, the nonlinear spring (modulus √ f − E 0 f E ; f is fibril stress, f is logarithmic fibril strain, and E is strain-dependent fibril network modulus) is connected in series with a linear dashpot (damping coefficient ) For a viscoelastic fibril i the Cauchy stress tensor is where z is the depth-dependent collagen fraction per solid volume, C is the ratio between primary and secondary fibrils, and ⊗ symbolizes dyadic product. In the collagen network there were 4 primary fibrils defining the Benninghoff arcadetype fibril formation and 13 secondary fibrils describing randomly oriented fibrils and collagen cross-links. e f is the current normalized fibril orientation vector The fluid flow within the tissue was incorporated according to Darcy's law where q is the flow flux, k is the current (hydraulic) permeability, and ∇p is the pressure gradient. The permeability was implemented as deformation dependent (van der Voet 1997) where k 0 is the initial permeability, e and e 0 are the current and initial void ratios, respectively, and the constant M describes how the permeability depends on the volumetric strain of the tissue.
At equilibrium conditions the Donnan osmotic swelling pressure gradient is where c F is the current depth-dependent FCD concentration, int , ext , ± int and ± ext are internal and external osmotic coefficients and internal and external activity coefficients, respectively (Huyghe et al. 2003), c ext is the external salt concentration (0.15 M), R is the molar gas constant (8.314 J/ mol K), and T is the absolute temperature (293.0 K). The stress due to chemical expansion is where a 0 and are material constants (Wilson et al. 2005b) and c − is the mobile anion concentration. The current depthdependent FCD concentration is modeled as a function of volumetric deformation where c F,0 is the initial depth-dependent FCD and n f,0 is the initial fluid volume fraction (i.e., porosity). Now, the total stress tensor in each material point is given by where totf is the sum of primary and secondary fibrils (4 + 13 = 17) and f is the chemical potential of water (Wilson et al. 2005a). The structural, compositional and material parameters of mature bovine articular cartilage are presented in Table 1.

Simulations and boundary conditions
Based on the cartilage lesion geometries found from clinical (Venäläinen et al. 2016) and in vitro (Orozco et al. 2018) studies and shown in Fig. 2a, a cylindrical 3D FE model with element type C3D8P representing an injured chondral plug (thickness h = 1 mm; radius r = 1.5 mm) was created along with an intact reference model (Fig. 2b). The FRPVES material model was implemented by using the user-defined material script (UMAT) subroutine in ABAQUS 6.14-2 (Dassault Systemes, Providence, RI, USA). Both models were subjected to repetitive 2-MPa axial ramp load in unconfined compression with a rigid impermeable compression platen. The duration of the ramp loading was 0.1 s, which together with the load magnitude represents typical values of the gait loading response (Gilbert et al. 2014;Kłodowski et al. 2016). Such loading values have been used in a previous computational study as well ). Unconfined compression was chosen due to its easy experimental setup and because it has been commonly used in other studies (Gratz et al. 2009;Szarko and Xia 2012;Li et al. 2013;Tanska et al. 2018). The lesion had the width of the explant since such lesions could be reproduced experimentally.
First, the cartilage was allowed to swell freely to reach mechanical equilibrium with the external environment (external salt concentration c ext = 0.15 M). The fluid was allowed to flow through the free edges during the free swelling step (pore pressure = 0): for the reference model the cylinder edges and the top surface, and for the injury model the lesion edges in addition to these. Subsequently, the compression platen was brought in contact with the surface (contact type: surface-to-surface; contact behavior: no separation and hard pressure-overclosure). The friction coefficient between the cartilage surface and the compression platen was set to 0.05 [cartilage-to-cartilage contact without meniscus (McCann et al. 2009)]. During these two  for comparison of computationally predicted FCD losses and in vitro findings by Orozco et al. (2018). c Outline of the proposed proteoglycan (PG) degeneration algorithm based on the decrease in the fixed charge density (FCD) of cartilage. While subjecting the explant in unconfined compression with an impermeable compression platen with 2-MPa uniaxial load exerted in 0.1 s, the FCD decrease was simulated based on either exceeded strains or stresses steps the bottom nodes were fixed to restrict translations of the tissue. Radial expansion of the tissue was allowed during the unconfined compression step. During the compression the fluid was allowed to flow only through the explant sides and lesion. After compressive loading, the mechanical analysis was conducted by characterizing the stress and strain measures followed by an evaluation of the FCD loss in the model (see Sect. 2.5). Next, a new FCD distribution was implemented in the model and a new analysis iteration was started using the same loading protocol. This process was repeated for 50 iterations (i.e., arbitrary time) to ensure the FCD distribution had reached an equilibrium state (Fig. 2c).

Mesh density
Simulations with denser mesh densities showed no substantial differences in maximum shear strain distributions (Fig. 3). The mesh with 4974 elements was used for the rest of the study since by doubling the mesh density the deviation of the peak strains was 0.7%. Also, the predicted distributions of decreased FCD contents with this mesh density were qualitatively and quantitatively the same (i.e., the same regions were degenerated with similar FCD loss) compared to denser meshes ( Fig. 3a-c). Additionally, the degeneration threshold and the degeneration rate (see Sect. 2.4) did not affect the mesh convergence (data not shown).

Proteoglycan degeneration
In the earlier studies including both clinical (Venäläinen et al. 2016) and experimental (Orozco et al. 2018) data, the found cartilage lesions (Fig. 2a) have been associated with changes in the mechanical behavior of the tissue. These changes caused by abnormal loading conditions and/or excessive mechanical stimuli may lead to chondrocyte death and PG loss (Fig. 2a, in vitro). Therefore, in the current study the PG degeneration and the subsequent FCD loss were based on the internal strain and stress distributions within the cartilage. Furthermore, the presence of lesions may modulate these strain/stress distributions locally, thus leading to altered FCD concentration levels.
In the current study, the FCD degeneration thresholds were obtained from the literature. The computationally predicted FCD losses were then compared to the in vitro findings by Orozco et al. (2018) (Fig. 2a). Furthermore, we assessed the sensitivity of the degeneration threshold, load magnitude and altered tissue properties (reduced collagen fibril integrity and increased permeability) on the simulated FCD loss.
The outline of the PG degeneration algorithm is shown in Fig. 2c. After obtaining the internal strain/stress distribution, the FCD decreases if the strain or stress measure exceeds the degeneration threshold ̂0 or ̂0 , respectively. Note that the strain-and stress-driven algorithms are separate, i.e., only one of them is used at a time. The strain-dependent degeneration rate ̂r and stress-dependent degeneration rate ̂r depend on the local strain and Based on this analysis, the mesh with 4974 elements was selected for computational efficiency since the strain distributions were mesh independent (peak strains deviated 0.7% compared to the mesh in a) and the resulting average local fixed charge density (FCD) distributions at different regions (see Fig. 2b) were qualitatively and quantitatively the same with the denser meshes (i.e., same regions were degenerated with similar FCD loss) stress level, respectively; the higher the strains/stresses, the faster the FCD concentration decreases via consecutive loading iterations of arbitrary time. Let ̂ and ̂ be the strain and stress values in an element (see Sect. 2.5 for more details), respectively, both calculated as the mean of strain/stress values at the integration points of the element. Then, according to an earlier study (Mononen et al. 2018), for the strain-driven algorithm the degeneration rate ̂r is and for the stress-driven algorithm Therefore, the amount of FCD for the iteration j in each element el in the strain-based algorithm is and in the stress-based algorithm The PG degeneration algorithm was implemented with a custom-made MATLAB script (R2017b, The Math-Works, Inc., Natick, MA, USA), which altered the FCD concentration in each element. The implemented initial FCD is affected by the volumetric deformation of the tissue during the mechanical analysis (Eq. 13). At the end of an unconfined compression step, the strains (which are affected by the deformation-modulated FCD) are imported into the degeneration algorithm along with the initial FCD concentrations of the current degeneration iteration. The FCD output from the degeneration algorithm (Eqs. 17 or 18) is then imported as a new FCD ( c F,0 in Eq. 13) for the next ABAQUS analysis iteration (50 iterations in total). Ultimately, the change of the FCD concentration has an impact on the Donnan osmotic swelling pressure gradient Δ and the chemical expansion stress T c (Eqs. 11 and 12). The decrease in the FCD has an effect on the total stress tensor tot (Eq. 14). For computational stability and avoiding convergence issues, the FCD concentration was not allowed to decrease to zero during the 50 degeneration iterations. Thus, the minimum allowed amount of the FCD concentration was set to 10% of the smallest initial FCD concentration of the model (found from the superficial layer).

Strain/stress measures and degeneration thresholds
Five strain measures and one stress measure were tested for the FCD depletion with the injury and intact reference model (Fig. 4a). In the PG degeneration algorithm, the degeneration thresholds ̂0 and ̂0 were chosen according to the literature. The total (integrated) strain and the Cauchy (true) stress tensors at the end of the ramp load were used for calculating the stress or strain measure in question. These tensors were obtained with UVARM and UEXTER-NALDB subroutines in ABAQUS. The following strain or stress measures were evaluated for PG degeneration: (i) Axial strain The absolute value of the axial direction (y)-component of the total strain tensor. The degeneration threshold for this non-invariant (dependent on the coordinate system) strain measure was ̂0 = 30% , which has been proposed to be injurious for the cartilage matrix ( . We have chosen to use "true octahedral shear strain" (Dixit and Dixit 2008), which has been used in other studies as well (Hosseini et al. 2014;Mononen et al. 2018). It can be written as where p,i are the principal strains of total strain tensor. The degeneration threshold was ̂0 = 17% as in Mononen et al. (2018). (iv) Maximum shear strain This invariant strain measure can be expressed as follows (Sibole and Erdemir 2012): where p,i are the principal strains of total strain tensor. The presence of defects has been shown to increase shear strain levels (Gratz et al. 2008;Venäläinen et al. 2016), which suggests that the maximum shear strain is a plausible parameter for PG degradation. The degeneration threshold was set to ̂0 = 30%, which has been linked with cartilage matrix damage ( i.e., the minimum principal strain value of total strain tensor. The presence of cartilage defects has been proposed to increase the minimum principal strains (Venäläinen et al. 2016). The degeneration threshold was chosen to be ̂0 = 30% (proposed threshold i.e., the maximum principal stress value of the Cauchy stress tensor. Here, the degeneration threshold is presumably associated with collagen network damage. The damaged collagen network does not retain PGs as well as a healthy network. Danso et al. (2014), Williamson et al. (2003) and Eleswarapu et al. (2011) have found experimentally that bovine cartilage ultimate tensile stress is 4.7 ± 2.1 MPa, 4.3 MPa and 7.0 ± 2.2 MPa, respectively. Therefore, the degeneration threshold was chosen as ̂0 = 5 MPa.
(22) = max p,1 , p,2 , p,3 , (a) PG degeneration algorithm (b) Further analysis based on strain/stress measure with altered conditions  Further analysis with altered conditions (i.e., different loading magnitudes, reduced collagen fibril integrity and increased tissue permeability) was conducted for the maximum shear strain-driven algorithm with a degeneration threshold consistent with earlier studies (D'Lima et al. 2001;Wilson et al. 2006;Venäläinen et al. 2016;Párraga Quiroga et al. 2017a). These analyses were carried out to test the sensitivity of the algorithm under loads ranging from a 0.5-MPa load to a 6-MPa load causing cell death (Fehrenbacher et al. 2003), and during simulated post-traumatic osteoarthritis conditions (halved initial ( E 0 ) and strain-dependent fibril moduli ( E ) and viscoelastic damping coefficient ( ))

Analysis of the in vitro and the simulated FCD loss
To compare the results of the degeneration algorithm to experimental findings, histological sections of two representative samples from the study by Orozco et al. (2018) are shown in Fig. 5a. In these experiments, sample sets were collected from patellofemoral grooves of 1-week-old calves.
The first set of the samples (n = 10) were injuriously loaded (50% strain at loading rate of 100%/s) to create defects prior to dynamic compressive loading (haversine waveform with 15% strain amplitude at 1 Hz for 1 h, followed by a 5-h rest period; repeat this loading cycle 4 times per day for 12 days; harvest of samples at day 12). Furthermore, the second set of the samples (n = 11) were injuriously loaded but not subjected to subsequent dynamic loading (harvest of The negative numbers show the average relative decrease in the FCD concentration at the cartilage appendage, tip and rim of the defect between the initial and simulated FCD concentrations samples at day 0). For this study, the relative FCD loss in the representative injured samples subjected to dynamic loading was calculated from the optical densities in relation to the representative injured sample without loading. In our computational models, the relative FCD losses were calculated locally at the cartilage appendage, tip and rim of the defect (see Fig. 2b) from one-element-thick slices halfway between the center and the edge of the sample with respect to the initial FCD distribution.

Effect of the degeneration threshold on the simulated FCD loss
Based on the degeneration threshold values from the literature, a few additional thresholds (higher and lower than literature thresholds) were chosen to determine the effect of the degeneration threshold on the resulting FCD distribution (Fig. 4a).
(i) Axial strain Besides the degeneration threshold of ̂0 = 30% , the other tested strain thresholds were 15% and 50%. These both thresholds were chosen according to a visual inspection of the strain distribution with ̂0 = 30%. In addition, the threshold value ̂0 = 15% was previously used by Mononen et al. (2018). (ii) Shear strain Thresholds of ̂0 = 30% (literature value), ̂0 = 15% and ̂0 = 50% were tested as in the previous case. (iii) Octahedral shear strain Besides the degeneration threshold of ̂0 = 17%, we also tested degeneration thresholds of ̂0 = 10%, ̂0 = 24% and ̂0 = 30%. The threshold of ̂0 = 10% was based on a visual examination of the strain distribution from the simulations with ̂0 = 17%, whereas the 30% of strain has been proposed to be injurious for chondrocytes (D'Lima et al. 2001) and to cause matrix damage (Wilson et al. 2006;Hosseini et al. 2014;Venäläinen et al. 2016;Párraga Quiroga et al. 2017a

Effect of the loading magnitude and altered tissue properties on the simulated FCD loss
In addition to investigating the effect of different strain and stress measures, we also tested how the loading magnitude and changed material properties affect the FCD loss (Fig. 4b). The maximum shear strain with the degeneration threshold ̂0 = 30% was chosen to drive the algorithm. First, in addition to 2-MPa load, four different loads (0.5 MPa, 1 MPa, 4 MPa and 6 MPa) were tested to see how the FCD distribution changes under different loading conditions. To justify the chosen loads, dynamic loading with 6 MPa may result in cell death in chondral explants (Fehrenbacher et al. 2003). Although these levels of stresses may be reached during everyday activities (Gilbert et al. 2014), they may cause alterations to matrix structure and impair the functioning of chondrocytes (Martin et al. 2004). Secondly, the effect of damaged collagen network and increased tissue permeability, both present in PTOA (Mäkelä et al. 2015), on the FCD loss was tested. To simulate the damaged collagen network we reduced the collagen fibril integrity by 50% [i.e., by decreasing the initial and straindependent fibril network moduli ( E 0 and E ) as well as the viscoelastic damping coefficient ( )]. For the increased tissue permeability, the initial permeability k 0 was increased by fivefold to characterize the effect of increased permeability on the FCD loss.

Comparison of the experimentally observed and the simulated FCD loss using different strain and stress measures
In the previous in vitro experiments (Orozco et al. 2018), optical density (i.e., an estimate of the FCD concentration) was low in the superficial zone and around the lesions. Furthermore, substantial cell death was observed near the defects. The FCD losses in the dynamically loaded injured samples compared to the injured sample without loading were approximately 50-60% near the cartilage appendage, tip and rim of the defect (Fig. 5a). The simulated qualitative local FCD concentrations near the lesion are shown in Fig. 5b-g for the different degeneration algorithms (the degeneration thresholds were obtained from the literature: axial strain ( ̂0 = 30% (D'Lima et al. 2001;Wilson et al. 2006;Hosseini et al. 2014;Venäläinen et al. 2016;Párraga Quiroga et al. 2017a)), xy-shear strain ( ̂0 = 30%), octahedral shear strain ̂0 = 17% (Mononen et al. 2018)), maximum shear strain ( ̂0 = 30%), minimum principal strain ( ̂0 = 30%) and maximum principal stress ( ̂0 = 5 MPa (Danso et al. 2014)). All strain measures except axial strain showed substantial local FCD depletion at the rim of the defect (approximately 50-60% decrease in the FCD content) and below the cartilage appendage at the tip of the defect (20-70%).
In more detail, for the non-invariant strain measures (i.e., dependent of the chosen coordinate system), the axial straindriven algorithm predicted local FCD loss only on the rim of the defect (45%), while the xy-shear strain-driven algorithm predicted 50% and 20% FCD loss at the rim and tip of the defect, respectively. For the invariant strain measures, the local FCD losses at the defect rim and tip were 49% and 20% for the octahedral shear strain, 59% and 69% for the maximum shear strain, and 48% and 36% for the minimum principal strain. The maximum principal stress-driven algorithm showed local FCD loss only at the rim of the defect (36%) and outside the assessed regions. None of the strainor stress-driven algorithms predicted FCD loss at the cartilage appendage.
The same differences between the tested strain and stress measures are also observed from the simulated depth-wise FCD profiles in Fig. 6. In the injury model, none of the strain measures showed depth-wise FCD loss at the cartilage appendage (Fig. 6, profile 1), below (profile 2) nor away from the defect (profile 4) compared to the depth-wise FCD profile of the intact reference model. Instead, all the strain measures predicted depth-wise FCD loss at the rim of the defect (profile 3). Below the cartilage appendage near the tip of the defect (profile 1), all strain measures except axial strain showed notable depth-wise FCD loss compared to the intact reference FCD profile.

Effect of the degeneration threshold on the simulated FCD loss
The decrease in the FCD concentration was sensitive to the chosen degeneration threshold. In general, decreasing the degeneration threshold increased the FCD losses and the borders of the regions of FCD loss. Increasing the degeneration threshold resulted in an unaltered FCD distribution in most of the cases. These results are shown in Fig. 7 for the axial strain, in Fig. 8 for the xy-shear strain, in Fig. 9 for the octahedral shear strain, in Fig. 10 for the maximum shear strain and in Fig. 11 for the minimum principal strain. With low degeneration thresholds, the strain-driven algorithms showed slight/nonexistent (xy-shear strain, octahedral shear strain) or substantial (axial strain, maximum shear strain, minimum principal strain) FCD loss near the explant rim in the injury and intact reference models. When the degeneration threshold for the maximum principal stress-driven algorithm was increased, the FCD concentration was decreased in smaller regions in the injury model and in the superficial zone of the center of the explant in the intact reference model (Fig. 12). More detailed explanation of the effect of the degeneration threshold on the predicted FCD losses can be found from the Supplementary Material.

Effect of the loading magnitude and altered tissue properties on the simulated FCD loss
The effect of different compressive load magnitudes was also tested (in addition to 2-MPa load) using the maximum shear strain algorithm with the degeneration threshold of ̂0 = 30% (Fig. 13). By decreasing the loading magnitude to 1 MPa, a negligible FCD loss occurred in a small area on the rim of the defect (Fig. 13b, iii, 1 MPa), as opposed to the 2-MPa case where a noticeable FCD loss was observed also near the tip of the crack in the transitional zone and upper deep zone (Fig. 13a, i, 2 MPa). When the load was further decreased to 0.5 MPa, no FCD loss occurred in the injury model (Fig. 13b, iii, 0.5 MPa). In the corresponding intact reference models, no FCD loss occurred (Fig. 13b, iv, 0.5 MPa and 1 MPa). When the compressive loading magnitudes were increased, the results were amplified compared to the 2-MPa case, i.e., the regions of FCD loss were larger and the FCD concentration was further decreased. The degeneration threshold ̂0 = 30% was also exceeded on the explant edges, subsequently leading to FCD loss in these areas as well (Fig. 13a, i, ii, 4 MPa and 6 MPa). Moreover, with 6-MPa load, the FCD concentration of the injury model was decreased in the deep zone near the explant edges and on the both sides away from the lesion (Fig. 13b, iii, 6 MPa).
We also investigated the effect of decreased integrity of the collagen fibrils and increased hydraulic permeability of the tissue on the FCD loss (Fig. 14). In the injury model, the maximum shear strains were higher for both of these cases, and thus, the areas of FCD loss were larger compared to the case with the reference conditions (initial fibril network modulus E 0 , strain-dependent fibril network modulus E , viscoelastic damping coefficient and initial permeability k 0 from Table 1). In both the intact reference model and the injury model the depth-wise FCD concentration was also decreased in the superficial zone near the explant rim compared to original conditions (Fig. 14b, iii and iv, reduced collagen fibril integrity and increased tissue permeability).

Summary
In this study we developed an algorithm for simulating an injury-related PG loss of articular cartilage by decreasing the local FCD concentration. Several potential mechanisms based on the local tissue strains or stresses were tested. The FCD loss was analyzed in a local (around the lesion) and depth-wise manner. We also investigated the sensitivity of the degeneration threshold, applied loading magnitude and altered material properties [reduced collagen fibril integrity    (Fig. 6). Nonetheless, all the invariant (i.e., independent of coordinate system) strain measures yielded similar FCD loss in the same regions around the lesion if the degeneration threshold was chosen optimally.

Strain-based degeneration algorithm(s)
In the strain-based algorithm (Fig. 4a, i-v), the FCD loss was driven by the internal strains of the tissue. It is known that high internal strains are related to chondrocyte apoptosis (D'Lima et al. 2001). This leads to locally decreased synthesis of new extracellular matrix constituents. However, this chondrocyte mechanotransduction process, which subsequently leads to local FCD loss, was not directly modeled in this study. Nevertheless, if the degeneration threshold ̂0 is selected based on the strain levels related to the cell apoptosis, we can relate the FCD loss predicted by the algorithm to the decreased synthesis of PGs. In addition, the algorithm also considers the failure of cartilage matrix, which has been proposed to occur at similar strain magnitudes ( ̂0 = 30%) (Wilson et al. 2006) as the chondrocyte apoptosis. The failure of cartilage matrix may potentially lead to collagen fibril damage, which was not directly modeled. Instead, if the degeneration threshold ̂0 is selected based on the strain (a) (c) (b) Fig. 7 a The axial strain distributions in (i) injury and (ii) intact reference models using the degeneration threshold of ̂0 = 30%. b The fixed charge density (FCD) distributions after 50 iterations for different degeneration thresholds ̂0 for (iii) injury and (iv) intact reference models. c Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (v) injury and (vi) intact reference models levels related to matrix failure, the FCD loss may be related to the decreased ability of the damaged collagen network to retain PG aggregates within the tissue.
All tested strain measures predicted no local or depthwise FCD loss at the cartilage appendage (Figs. 5, 6, profile 1) as opposed to in vitro findings (Fig. 5a). The reasons for this could lie in folding of the appendage into the lesion instead of deforming, the fact that crack propagation was not modeled in the current study, different loading protocols and different material properties (discussed more in Sect. 4.7). On the contrary, at the tip of the defect below the cartilage appendage, the local and depth-wise FCD contents were decreased substantially with all strain measures except the axial strain (Figs. 5 and 6, profile 1). The FCD loss near the tip of the defect might be due to the folding of the cartilage appendage into the lesion, while the internal axial strains do not increase substantially. The largest decrease in the local (approximately 70% decrease in the FCD) and depthwise FCD contents occurred below the cartilage appendage using the maximum shear strain-based algorithm (Figs. 5e, 6d, profile 1). This strain measure exhibited the depth-wise FCD loss on the side of the cartilage appendage (profile 1) at the normalized depths from 0.3 to 0.4 (~ in the deepest 1/3 of the lesion). However, with the minimum principal strain-based algorithm, the local FCD loss near the tip of the defect was 36% (Fig. 5f) and in profile 1 the depth-wise FCD loss was located only at the normalized depth of 0.4 (Fig. 6e). This indicates that the maximum shear strain reached injurious levels in greater number of elements than the minimum principal strain in profile 1. Furthermore, the (a)  Fig. 8 a The xy-shear strain distributions in (i) injury and (ii) intact reference models using the degeneration threshold of ̂0 = 30%. b The fixed charge density (FCD) distributions after 50 iterations for different degeneration thresholds ̂0 for (iii) injury and (iv) intact reference models. c Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (v) injury and (vi) intact reference models minimum principal strain-based algorithm predicted smaller depth-wise FCD content in profile 1 (at the normalized depth of 0.3) compared to the axial strain-, the shear strain-and the octahedral shear strain-based algorithms (Fig. 6a, b, c and e). On the rim of the defect, the local decrease in the FCD was the largest with the maximum shear strain-based algorithm (approximately 60% decrease in the FCD, Fig. 5e) compared to the other strain-based algorithms (around 50% decrease, Fig. 5). The depth-wise FCD profiles on the rim of the lesion (Fig. 6, profile 3) were similar between all tested strain measures, and the profiles also resemble each other below the lesion (profile 2) and in the intact area away from the lesion (profile 4).
The octahedral shear strain (invariant strain measure) resulted in similar quantitative local FCD losses (approximately 20% decrease in the FCD concentration compared to the initial FCD concentration at the tip of the lesion, Fig. 5c and d) and depth-wise FCD losses as the shear  Fig. 6b, c). However, the octahedral shear strain (and thus, the related degeneration threshold) is more difficult to determine experimentally compared to the "regular" shear strain. Thus, the degeneration threshold for octahedral shear strain ( ̂0 = 17%) was obtained from another computational study (Mononen et al. 2018). For instance, the degeneration threshold for the "regular" shear strain ( ̂0 = 30%) used in this study is related to chondrocyte apoptosis (D'Lima et al. 2001). However, the fact that the shear strain was excessive is strongly related to the chosen geometry and coordinate system. Thus, it could be more advantageous to use invariant strain measures which do not depend on the chosen coordinate system, such as the maximum shear strain or the octahedral shear strain.
The PG loss has been observed to occur focally around lesions (Moeini et al. 2014;Mohanraj et al. 2014;Thomas et al. 2017;Orozco et al. 2018). Consistent with these findings, the largest decrease in the FCD concentration was observed around both sides of the lesion with the maximum shear strain-and minimum principal strain-based algorithms (Figs. 5e, f vs. 6d, e). Compared to the study by Orozco et al. (2018), the estimated experimental FCD losses were around 50-60% at the tip and rim of the defect (Fig. 5a), whereas the FCD losses predicted by the maximum shear strain-based algorithm were around 60-70%. The minimum principal strain-based algorithm predicted slightly smaller FCD losses (35-50%) at the tip and rim of the defect (Fig. 5e, f). For both of these strain-based algorithms, the degeneration was The fixed charge density (FCD) distributions after 50 iterations for different degeneration thresholds ̂0 for (iii) injury and (iv) intact ref-erence models. c Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (v) injury and (vi) intact reference models associated with the experimentally approximated degeneration threshold ( ̂0 = 30%) (D'Lima et al. 2001). Altogether, with this degeneration threshold, the most comprehensive FCD loss was observed with the maximum shear strain-and minimum principal strain-based algorithms, which coincides with the earlier observation of the elevated shear and minimum principal strains near cartilage defects (Venäläinen et al. 2016). These strain measures are strain invariants (i.e., independent of the coordinate system) and are able to replicate the FCD loss around the lesion, which is also observed experimentally (Moeini et al. 2014;Mohanraj et al. 2014;Thomas et al. 2017;Orozco et al. 2018). Thus, they could be potential strain measures for helping clinicians to decide which defects possess a high risk (due to its location and size (Peña et al. 2007;Venäläinen et al. 2016)) to cause tissue degradation. In addition, the octahedral shear strain-based algorithm exhibited qualitatively similar FCD loss around the lesion (in the same regions, i.e., on the defect rim and at the tip of the defect; see Figs. 5d, 6c), but its degeneration threshold is more challenging to determine directly from experiments. The results of the xy-shear strain were satisfactory only for this unconfined compression geometry. However, as xy-shear strain is dependent on the chosen coordinate system (which is important factor if the sample is rotating) and geometry, it is probably too simple for the complex in vivo loading conditions. Instead, the maximum shear strain is an invariant alternative to the xy-shear strain. The axial strain produced notable FCD loss only on the rim of the lesion (Figs. 5b, 6a). Decrease in FCD concentration % Fig. 11 a The minimum principal strain distributions in (i) injury and (ii) intact reference models using the degeneration threshold of ̂0 = 30%. b The fixed charge density (FCD) distributions after 50 iterations for different degeneration thresholds ̂0 for (iii) injury and (iv) intact reference models. c Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (v) injury and (vi) intact reference models Taken together, when the degeneration threshold was selected based on strain-related chondrocyte death and matrix failure, the maximum shear strain produced the largest FCD losses which are similar compared to the experimental FCD losses (Orozco et al. 2018). This conclusion is supported qualitatively by the comparable regions of the FCD loss (low optical density) (Fig. 5) and quantitatively by the most extensive local FCD losses near the lesion compared to the experimental findings (Fig. 5) and the other tested strain measures (Figs. 5, 6). However, our algorithm did not predict the FCD loss in all locations around the lesion, for example at the cartilage appendage (discussed in Sect. 4.7). We underline that the FCD loss would be similar for all the tested invariant strain measures (i.e., the octahedral shear strain, the maximum shear strain and the minimum principal strain) if the degeneration threshold could be chosen individually for each strain measure such that the FCD loss would be localized to the neighborhood of the defect. Naturally, these different degeneration thresholds should be determined experimentally. Furthermore, the degeneration threshold most likely depends on the status and initial material properties of the cartilage. For example, Orozco et al. (2018) used maximum shear strain-based algorithm with the degeneration threshold of 50% to match the computational model with experimental findings in calves. In the current study, we used material parameters representing mature bovine cartilage and the threshold used by Orozco et al. did not lead to similar FCD losses compared to the in vitro experiments (Fig. 10c, v, ̂0 = 50%). However, this is not the only explanation for the differences between the current model and the in vitro findings (see Sect. 4.7). Thus, the degeneration threshold for human tissue is most likely different. The value of the degeneration threshold and the degeneration parameter (for example, failure strain, stress, fluid flow) might also differ for different tissue degradation mechanisms, such as cell death, flow of fragmented PGs out of  Fig. 12 In the stress-based algorithm, excessive levels of maximum principal stress presumably damage the collagen network, thus impairing its ability to retain proteoglycans within cartilage. a The maximum principal stress distributions in (i) injury and (ii) intact reference models using the degeneration threshold of ̂0 = 5 MPa. b The fixed charge density (FCD) distributions after 50 iterations for different degeneration thresholds ̂0 for (iii) injury and (iv) intact reference models. c Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (v) injury and (vi) intact reference models the tissue, collagen failure. For instance, cartilage damage and degeneration might be directly related to the ultimate tensile strain of individual collagen fibrils (Hosseini et al. 2014), which could be considered to remain constant. Furthermore, the degeneration threshold may vary between different subjects due to aging, similarly as the ultimate tensile failure stress (Kempson 1982). These are naturally important aspects for subject-specific modeling.

Stress-based degeneration algorithm
The stress-based algorithm (Fig. 4a, vi) considered the FCD loss via the excessive maximum principal stresses, which presumably impair the ability of the collagen network to retain the PGs within cartilage.
On the contrary to experimental findings, there was no local FCD loss in and below the cartilage appendage (Fig. 5g). This was also seen from the depth-wise FCD profiles (Fig. 6f). The stresses elevated near the rim of the defect as in the previous computational study (Venäläinen et al. 2016), but they also increased at the intact adjacent tissue near the lesion (Fig. 12a, i, ii). As opposed to the  Fig. 13 a The fixed charge density (FCD) distributions after 50 iterations in (i) injury and (ii) intact reference models using different load magnitudes. The maximum shear strain with the degeneration thresh-old of ̂0 = 30% was used for the FCD depletion. b Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (iii) injury and (iv) intact reference models strain-based algorithm with invariant strain measures, the simulation results of the stress-based algorithm did not coincide with the hypothesis where the FCD losses would presumably occur only around the defect. The reason for the elevated stresses in large regions could lie in the fluid pressurization due to high loading rate. Altogether, although the maximum principal stress is suggested to be able to simulate collagen degeneration Liukkonen et al. 2017), it might not be a feasible mechanism for PG degradation.

Effect of the degeneration threshold on the simulated FCD loss
The results of the effect of the degeneration threshold on the resulting FCD distribution are depicted more in depth in the Supplementary Material. For the axial strain-, xyshear strain-, maximum shear strain-and minimum principal strain-based algorithms, the threshold of ̂0 = 30% was presumably appropriate, since the damage was constrained to the vicinity of the lesion. When the threshold was decreased to ̂0 = 15%, the areas of FCD loss became larger, but were still constrained to the adjacent tissue of the lesion with xy-shear strain. For example, when the lower degeneration threshold was used for the maximum shear strain in the intact reference model, the FCD concentration decreased on the explant rim and beneath the superficial zone near the explant edges ( Fig. 9c, vi, ̂0 = 15%). We are not aware of any experimental study showing similar FCD loss. Together, these findings suggest that the degeneration threshold of ̂0 = 15% may be too low for the maximum shear strain to predict the strain-related FCD loss. In the literature, a lower shear strain threshold has been proposed for cell death ( ̂0 = 7-12%) (Bonnevie et al. 2018) and presumed subsequent FCD loss. However, in the view of the current study, with the axial strain, maximum shear strain and minimum principal strain the other areas besides the near vicinity of the lesion were also damaged. This could indicate that the degeneration threshold value below ̂0 = 15% is too low for these strain measures. For the octahedral shear strain, the degeneration threshold of ̂0 = 17% decreased the FCD concentration on the both sides of the lesion (Fig. 7b), which did not occur for the two larger thresholds ( ̂0 = 24% or ̂0 = 30%). The threshold of ̂0 = 10% is most likely too  Fig. 14 a The fixed charge density (FCD) distributions after 50 iterations in (i) injury and (ii) intact reference models for the cases with either the reduced collagen fibril integrity or increased tissue permeability. For the case with reduced collagen fibril integrity, the viscoelastic damping coefficient , initial fibril modulus E 0 and straindependent fibril modulus E were reduced by 50%. For the case with the increased initial permeability ( k 0 ), k 0 was five times higher compared to the case without any alterations. The maximum shear strain with degeneration threshold of ̂0 = 30% was used for the FCD depletion. b Corresponding relative decrease in the FCD concentrations compared to the initial concentrations for (iii) injury and (iv) intact reference models small, as we observed minor FCD loss on the edges of the intact reference and injury models. This FCD loss was seen in the experiments, but was supposedly caused by sample preparation and fluid flow (Orozco et al. 2018). The optimal threshold could be between ̂0 = 10% and ̂0 = 17%; possible values could be ̂0 = 15% proposed in Mononen et al. (2018) or ̂0 ≈ 12% in Bonnevie et al. (2018), and it is possibly modulated by the age of the tissue.
For the stress-based degeneration algorithm, the degeneration threshold of ̂0 = 5 MPa was probably too small since the FCD decreased in large areas in both the injury and intact reference models primarily in the superficial zone (Fig. 12b, iii and iv, ̂0 = 5 MPa) in contrast to the experimental (Moeini et al. 2014;Mohanraj et al. 2014;Orozco et al. 2018) and computational (Hosseini et al. 2014;Mononen et al. 2018) studies, where the PG degradation of the injured tissue occurred around the lesion in both the superficial and transitional zones. However, the threshold of ̂0 = 7 MPa, as suggested in the literature Liukkonen et al. 2017), yielded results which conformed better with the experiments by Orozco et al. (2018) in a sense that the damage was concentrated to the defect rim and nearby regions in the injury model, whereas the intact reference model was damaged only from a minor area.

Effect of the loading magnitude and altered tissue properties on the simulated FCD loss
When the cartilage explant was subjected to different loads, the FCD loss increased as the magnitude of the load increased (Fig. 13). This load-dependent FCD loss is consistent with an earlier experimental study where the amount of dead cells increased with peak compressive load (Fehrenbacher et al. 2003), which may lead to a decrease in the local FCD due to decreased synthesis. Furthermore, when the two commonly observed changes in PTOA [i.e., reduced collagen fibril integrity (Punzi et al. 2016) and increased tissue permeability (Sauerland et al. 2003)] were tested, the FCD loss was increased in both cases, as expected. Similar changes in tissue properties and decrease in GAG coverage (associated with FCD loss) have been observed in experimental models of PTOA (Fischenich et al. 2017). In the injury model, the increased permeability led to more FCD loss near the lesion than decreased integrity of the collagen fibrils. On the other hand, in the intact reference model at the lateral edges of the explant, halving the collagen fibril integrity decreased the FCD concentration more than increasing cartilage permeability (Fig. 14b, iv, reduced fibril integrity and increased tissue permeability). Both changes in tissue properties ultimately reduce the load-bearing capability of cartilage. These changes result in higher strains in the tissue and subsequently lead to greater FCD loss.

Limitations
In the simulations there was no relaxation period, since the algorithm decreased the FCD concentration based on the local strains at the end of the ramp load (highest load at 0.1 s). This loading protocol and timing approximates the in vivo conditions before and during the first peak load of the stance phase of gait; especially at the peak (~ 0.1 s) the strains may exceed the degeneration threshold (Venäläinen et al. 2016) when assuming that the stance phase lasts ~ 0.8-1 s. However, if the FCD was decreased based on the local strains after a relaxation period, the strains would have increased slightly (data not shown) compared to the strains at the end of the ramp load. After relaxation there is less fluid in tissue to carry the load, thus resulting in increased strain levels which would have led to greater predicted FCD losses.
The degeneration rate ̂r should be optimized to some timescale according to experimental results instead of using the simple equation as in the current study, so that the model could be linked to physical timescale instead of arbitrary time. For example, one iteration could be scaled to correspond the FCD loss of one loading step in an in vitro experiment. This would bring the algorithm a step closer to in vivo use; then, the degeneration rate ̂r could be scaled, for example, to the FCD loss of half a year of walking.
One thing to consider in the model would be to initially increase the fluid fraction near the lesion, since injurious loading has been suggested to lead to tissue swelling (Torzilli et al. 1999). Increased tissue hydration has also been shown to lead to increased permeability (Gu and Yao 2003). This would ultimately affect the local osmotic pressure near the crack and change the total stress and strain tensors, thus probably yielding a different FCD distribution. However, it is unclear how fast the initial fluid fraction increases after lesion formation. Thus, we did not test this.

Other factors affecting the FCD loss and comparison to literature
Other factors affecting the FCD depletion could be increased tissue permeability and fluid velocity near the defects. The latter approach has some computational support by Orozco et al. (2018), but experimental studies on the effect of fluid velocity should be carried out. However, it could be possible that the effects of fluid velocity and the effects of the maximum shear strain from the current study could be combined in order to obtain a more comprehensive model for the FCD loss.
The interaction between tissue hydration, permeability and FCD affects profoundly the intra-tissue strain distributions. This interaction has also a major influence on tissue health and depth-wise OA progression. If no electrical current flows through the tissue, such as in cartilage, the open-circuit permeability is regulated by at least two mechanisms: (1) the mechanical effect where the porosity of cartilage changes as a function of tissue deformation or alteration in the GAG content (i.e., intrinsic closed-circuit FCD-independent permeability) and (2) the FCD-related electro-kinetic coupling effect in which the fluid flow is affected by electrical potential differences (electro-osmosis phenomenon) (Gu and Yao 2003). The open-circuit permeability is more sensitive to changes in water content than changes in the FCD (Gu and Yao 2003). Tissue hydration is in turn related to the applied load/strain. External mechanical loading decreases tissue hydration, which decreases porosity and thus permeability. Hence, the permeability is strain dependent (van der Voet 1997; Gu and Yao 2003). High water content is related to high permeability which leads to higher rate of deformation (Gu and Yao 2003) and higher strains. In early OA, especially in the superficial zone, both high water content and elevated strains are found along with disorganized collagen fibrils and decreased FCD content (Saarakkala et al. 2010). On the other hand, high FCD content is related to small permeability (Maroudas 1975;Gu and Yao 2003) and water flux (Gu and Yao 2003) and, thus, to smaller strains. Sufficient hydration and fluid pressurization during loading are also important, which is supported by the finding that dehydrated cartilage has an increased risk of rupturing compared to well-hydrated cartilage (Fick and Espino 2011). Furthermore, high FCD concentration ensures more efficient rehydration of the tissue after loading . Taken together, both high FCD concentration and normal tissue hydration are essential to keep the intratissue strains at a nondestructive level. This behavior is also reflected in our model, since low-FCD regions experienced larger strains compared to high-FCD regions (due to the lesion but also during the later degeneration iterations due to decreased FCD content). In our model, the superficial zone was only slightly deformed which is accounted for by sufficient hydration, high fluid pressurization due to collagen fibril orientations (Pierce et al. 2016) and their viscoelastic nature (Párraga Quiroga et al. 2017b).
In our model, the strains elevated primarily near the lesion for all tested strain measures (Figs. 7,8,9,10,11a,12a). Similar increase in strains near the lesions was also seen by Gratz et al. (2009) andZevenbergen et al. (2018) in pairs of osteochondral explants and by Venäläinen et al. (2016) in the knee joint. These studies support the strain distribution findings in our study. Venäläinen et al. also found that all strains except axial strains increase beneath the lesion. We observed similar behavior of axial strains, but in our model the other strain measures were not elevated beneath the lesion. This might be due to the more complex loading and geometry (gait cycle and knee joint) as well as lesion shape used by Venäläinen et al. In contrast to our findings, in an unconfined compression experiment by Moo et al. (2014), tissues far away from the cut (full-thickness cuts were created by cutting half the tissue away perpendicular to the tissue surface) experienced larger deformation (axial strain in superficial zone) than tissues close to the cut. This was seen experimentally by tracking cells and computationally with an FE model. Possible explanations for this difference are the near-equilibrium conditions (final displacement was held for 20 min) and different lesion geometry used by Moo et al. (full-thickness cut). On the contrary to our model, the cut might have caused depth-wise PG loss at the sample edge, which might have led to larger deformation in the deep zone and subsequently to smaller deformation in the superficial zone on the edge of the cut.
The FCD losses were similar to the experimental findings by Orozco et al. (2018) (Fig. 5a), but as opposed to these findings, our model did not predict any FCD loss near the cartilage appendage. Orozco et al. did not either predict FCD loss at the cartilage appendage with maximum shear strain-driven algorithm. The reasons for the FCD loss differences in the current study and the in vitro findings are different lesion geometries, different material properties of the samples [young calve (Orozco et al. 2018) vs. mature bovine], different loading protocols and, especially in our model, folding of the appendage into the lesion instead of deforming considerably because of fluid pressurization due to collagen fibril orientations (Pierce et al. 2016). Hosseini et al. (2014) and Párraga Quiroga et al. (2017a) used a damage progression model based on mechanical loading, which included the separate degradation of both collagen fibril network and PG-rich ground substance/nonfibrillar matrix (by decreasing the corresponding network or substance moduli). The driving factor for cartilage degeneration was the deviatoric strain (octahedral shear strain) with the non-fibrillar matrix damage initiation threshold set to 30% of strain. Indentation of the tissue resulted in PG damage in the superficial and transitional zones only near the indenter. This is comparable to the strain-based algorithm when the invariant strain measures (for example, the octahedral shear strain) are used, so the strain-based algorithm in the current study provides a great alternative to the damage model of the non-fibrillar matrix by Hosseini et al. and Párraga Quiroga et al. The joint level simulations by Mononen et al. (2018) also included the collagen network degeneration (by reducing the fibril network moduli) as well as the PG depletion (by reducing the non-fibrillar matrix modulus). The PG depletion was also coupled with alterations in cartilage permeability when the joint was under a physiologically relevant gait loading. Similarly, they used the deviatoric strain-based degeneration algorithm for the non-fibrillar matrix (the degeneration threshold of 17%) and found that the non-fibrillar matrix modulus (~ estimate of PG content) was decreased at the superficial and transitional zones. Again, qualitatively similar results were found near the lesion with the strain-based algorithm driven by the invariant strain measures.
The algorithms investigated in this study were only concentrated on the FCD loss due to a mechanical stimulus and excluded the effect of any biological processes such as inflammation or altered PG synthesis due to injurious loading of the chondrocytes. The presence of biochemical stimuli is evident in PTOA, and the damage models by Kapitanov et al. (2016) and Gardiner et al. (2016) included them. The former model combined cellular (healthy, diseased and dead cells), chemical (pro-and anti-inflammatory cytokines) and mechanical (tissue internal strains which induce the biochemical processes) components. However, the model excluded the regulation of matrix constituent synthesis rate by the state of cells. Furthermore, the study used a simple external pressure load that affected only the initial conditions. The latter model implemented the cell viability and related PG and collagen synthesis. The primary factor driving the amount of PGs is the ratio of ECM synthesis and degeneration in the equilibrium state instead of the direct mechanical damage. While these models are very sophisticated, neither one implemented the heterogeneous and anisotropic structure of the articular cartilage. By combining a material model including this intrinsic cartilage structure with similar biochemical mechanisms, the internal strains would change. This would affect depth-wise cell death and cartilage degeneration rate. For example, the maximum shear strain-driven algorithm could be incorporated into that model. If successful, such a model could encompass the complex biomechanical and biochemical environments and estimate PTOA progression.

Future
The proposed algorithm requires a thorough in vitro validation in the future. On top of the biomechanical stimuli that this algorithm considers, the biological repair responses and biochemical factors such as aggrecan synthesis, inflammatory responses (e.g., due to cytokines TNFα, IL-1α) (Kar et al. 2016) and shear-induced release of reactive oxygen species (Martin et al. 2004) should be taken into account in the future when aiming toward more accurate and validated in vitro and later in vivo models. Furthermore, the changes in tissue permeability due to the PG loss, collagen degeneration (as implemented in (Hosseini et al. 2014;Kar et al. 2016;Mononen et al. 2018)) and reorientation (as implemented in ) could be included to obtain even better estimation of the cartilage degeneration processes on a long-term basis.
The current model did not include the possible changes in subchondral bone properties due to PTOA, such as an altered density and plate thickness (Myller et al. 2017). These may affect the internal strain/stress distributions, especially near the cartilage-bone interface and calcified cartilage. According to the developed algorithm no FCD loss was observed in the deep cartilage directly below the lesion, which might be different in vivo (Myller et al. 2017). In reality, the subjectspecific in vivo loading conditions are more complex and they are not fully simulated in our simplified unconfined and axial compression schemes. Thus, more complex loading schemes including shear should be investigated in the future.

Conclusions
In conclusion, the strain-based FCD degeneration algorithm driven with invariant strain measures was able to predict the FCD distribution in the damaged articular cartilage tissue. Consistent with our hypothesis and also observed in the literature, the algorithm predicted the FCD loss to occur in the vicinity of the lesion. Furthermore, if the degeneration threshold was selected based on the strain levels related to cartilage matrix damage and chondrocyte apoptosis, the most comprehensive FCD losses around the lesion were obtained with the maximum shear strain. These local FCD losses were in accordance with the current experimental knowledge. However, the strain-based algorithm needs extensive in vitro validation in the future. The mechanical stimulus related to the integrity of the PG matrix is a key factor for the adaptation of cartilage during PTOA. On top of the biomechanical aspects, also the biochemical factors (such as PG synthesis and effect of pro-inflammatory mediators) should be considered in more refined PTOA progression models. Nevertheless, the strainbased FCD degeneration algorithm could be a part of a larger modeling framework to predict progression of PTOA and identify high-risk defects. The algorithm could help in decision making for early clinical interventions and potentially open up new pathways for early clinical applications.