## Abstract

The osmotic pressure in articular cartilage serves an important mechanical function in healthy tissue. Its magnitude is thought to play a role in advancing osteoarthritis. The aims of this study were to: (1) isolate and quantify the magnitude of cartilage swelling pressure in situ; and (2) identify the effect of salt concentration on material parameters. Confined compression stress-relaxation testing was performed on 18 immature bovine and six mature human cartilage samples in solutions of varying osmolarities. Direct measurements of osmotic pressure revealed nonideal and concentration-dependent osmotic behavior, with magnitudes approximately 1/3 those predicted by ideal Donnan law. A modified Donnan constitutive behavior was able to capture the aggregate behavior of all samples with a single adjustable parameter. Results of curve-fitting transient stress-relaxation data with triphasic theory in febio demonstrated concentration-dependent material properties. The aggregate modulus H_{A} increased threefold as the external concentration decreased from hypertonic 2 M to hypotonic 0.001 M NaCl (bovine: $HA=0.420\xb10.109\u2009MPa$ to $1.266\xb10.438\u2009MPa$; human: $HA=0.499\xb10.208\u2009MPa$ to $1.597\xb10.455\u2009MPa$), within a triphasic theory inclusive of osmotic effects. This study provides a novel and simple analytical model for cartilage osmotic pressure which may be used in computational simulations, validated with direct in situ measurements. A key finding is the simultaneous existence of Donnan osmotic and Poisson–Boltzmann electrostatic interactions within cartilage.

## 1 Introduction

Despite the critical role of osmotic pressure in biological tissue mechanics, the precise magnitude and concentration dependence of the osmotic pressure within osmotically active tissues remains elusive. In articular cartilage, the negatively charged glycosaminoglycan (GAG) chains enmeshed within the collagen network of the extracellular matrix (ECM) give rise to the well-known fixed charge density (FCD) of cartilage. When immersed in an electrolyte solution, such as the interstitial fluid which fills the pore space of cartilage, the abundance of fixed charges produces a swelling pressure within the tissue. ECM homeostasis is a balance between the osmotic swelling pressure exerted by the GAGs and tensile stresses developed in the collagen network [1], therefore producing a tissue which is prestressed even in situ.

Early osteoarthritis (OA) produces damage to the collagen network of cartilage without a significant change in the GAG content [2], disrupting this homeostasis. Cartilage then swells under the influence of the insufficiently restrained GAGs, producing a locally softer tissue with increased water content [2]. These locally softer regions of tissue are further susceptible to damage during physiological loading, and these sites are often loci where degeneration progresses. The cartilage swelling resulting from this cascade is clinically recognized as an early sign of OA [2,3]; this swelling has been observed in recent cartilage fatigue and damage studies [4,5] which sought to study the mechanically mediated damage process in cartilage. In the study by Durney et al., physiologic sliding conditions produced formation of large swollen blister-like sacs on the cartilage surface within the loaded region, with further loading leading to complete delamination of the superficial zone remarkably reminiscent of arthritic human cartilage [5]. Several recent studies underscore the importance of understanding and predicting the processes leading to damage and OA progression in cartilage [6–9]. A more complete understanding of such complex interactions in ECM mechanics is crucial in understanding mechanically mediated diseases such as OA.

Osmotic pressure is the pressure produced in a fluid compartment due to an imbalance in solute concentrations with an adjoining fluid compartment. Donnan osmotic pressure is the osmotic pressure produced when this imbalance is caused by electrostatic charges. In articular cartilage, the electrostatic charges are carried by ions in the interstitial fluid and proteoglycans enmeshed within the porous solid matrix. Theoretical models of cartilage mechanics have long recognized the central importance of osmotic and electrostatic forces in mediating the mechanical response of cartilage. The classical triphasic theory of cartilage [10] uses continuum mixture theory to model cartilage as a mixture of a charged solid phase, neutral fluid phase, and monovalent salt ions such as Na^{+} and Cl^{−}. The osmotic interactions enter this framework as a Donnan swelling pressure [11], caused by an imbalance of anion and cation concentrations in the tissue versus its surrounding bath [12]. In contrast, Grodzinsky et al. [13–16] have advocated for a hybrid approach termed continuum electromechanics, which combines electromagnetism and continuum mechanics to produce theories of the Poisson–Boltzmann type, where swelling is caused by charge-to-charge repulsion along parallel GAG side chains [17]. Microstructural constitutive models have been proposed to connect the two phenomena, by deriving the Donnan model from electrostatic interactions in appropriate limiting cases [18]. Our study works within the framework of triphasic theory and considers osmotic pressure to be of the Donnan type: The osmotic pressure arises from an osmolarity imbalance between the tissue interstitial fluid and the external bathing environment, caused by attractive and repulsive electrostatic forces between proteoglycan fixed charges and mobile ions.

Donnan law relies on constitutive models for the chemical potential of water (the solvent) and ions and is strongly dependent on empirical coefficients describing deviations from chemically ideal conditions. Due to nonlinear physicochemical effects, the Donnan osmotic swelling pressure in cartilage cannot be predicted accurately solely from the tissue's FCD using ideal Donnan law. Therefore, a direct experimental characterization of the Donnan pressure inside articular cartilage in situ is an essential prerequisite to a validated theory of cartilage damage. Despite extensive experimental [19–25] and theoretical [12,15,26–28] studies, the precise magnitude of Donnan osmotic pressure, and its relationship to cartilage composition, remains unknown. Studies on solutions of extracted GAGs have demonstrated that they behave in a nonideal manner [19,21,29,30], with ideal Donnan behavior estimated to over-predict swelling pressure by up to a factor of two [29]. However, it appears measurements of GAGs in solution do not accurately predict the Donnan pressure of cartilage in situ [31–33]; further, due to the general nonequivalence between mechanical and osmotic loading [32], such results must be interpreted carefully. To date, no studies have provided a validated measurement and constitutive model of this pressure within articular cartilage. Consequently, most theoretical and computational models assume ideal conditions (e.g., Refs. [24] and 34]).

Another open question concerns the effect of ion concentrations on material properties of the porous cartilage matrix, including its compressive modulus and permeability to interstitial fluid flow. It has been known for decades that cartilage stiffens significantly in hypotonic solutions and softens in hypertonic solutions (e.g., Ref. [23]). Theoretical correspondence principles have been proposed to include the effect of external salt concentration on tissue moduli [24,27], but currently lack experimental validation. In one of the most detailed experimental studies employing triphasic theory, Lu et al. [35] performed creep indentation on cartilage samples immersed in isotonic 0.15 M NaCl and hypertonic 2 M NaCl. Using triphasic theory and ideal Donnan behavior, these authors were able to find a linear relationship between the theoretically predicted FCD and the measured FCD of tested samples and demonstrated that osmotic effects influence the measured material properties. They concluded that there exists a set of intrinsic material properties, i.e., those measured in hypertonic conditions, which may be related to apparent properties measured at other concentrations through the osmotic pressure. The work of these authors, though comprehensive, only studied two external bathing concentrations and reported a low average FCD of 149 ± 59 mEq/L; these results may have precluded significant nonideality and might not generalize to higher FCDs present in bovine tissue [36], which has been recently used for cartilage damage studies (e.g., Refs. [5] and [9]). Additionally, the theoretical FCD overpredicted the measured values by more than a factor of two, indicating that some of the theoretical assumptions may not have been met [35].

To address these open questions, this study aimed to directly measure and characterize the concentration- and strain-dependent osmotic pressure and material properties in articular cartilage. To achieve these aims, this study (1) isolated and quantified the magnitude of cartilage swelling pressure in situ through a comprehensive experimental protocol for confined compression stress-relaxation and (2) identified the effect of ion concentration on so-called intrinsic material parameters using a triphasic framework that explicitly accounts for osmotic behavior. These experiments allowed direct measurement of the osmotic pressure of cartilage in situ and correlation with GAG content in the tissue. To examine nonideal behavior, bovine deep zone cartilage was used, thus ensuring a high FCD relative to other similar experiments and promoting nonideal conditions. An additional study on human cadaveric tissue examined the applicability of bovine results to mature human tissue. The experimental design measured the osmotic pressure on the same sample at six different bath concentrations, ranging from hypertonic 2 M to hypotonic 0.001 M, using a strictly defined reference frame in a procedure adapted from Eisenberg and Grodzinsky [23]. For stress-relaxation responses, two prescribed strain values were investigated alongside the six concentrations, thus providing a broad range of applicability. The results can inform computational models and experimental designs within the cartilage and intervertebral disc community, and further our aims of investigating the mechanically mediated damage process producing OA in articular cartilage.

## 2 Materials and Methods

### 2.1 Mathematical Modeling

#### 2.1.1 Confined Compression.

*σ*in triphasic theory reduces to

*p*is the fluid pressure inclusive of osmotic effects [34] and

*σ*is the elastic stress; the constitutive model for

^{e}*σ*depends on strain (e.g., the axial normal engineering strain

^{e}*ε*) and the concentrations of NaCl salt ions.

^{1}For this 1D analysis, based on the quasi-static mixture momentum balance, the mixture stress

*σ*is uniform along the axial (compression) direction. In the steady-state response to loading, the fluid pressure subsides to its osmotic contribution, $p=\pi $, and the concentration of salt ions in the tissue may be expressed as a function of the tissue strain

*ε*and the external bath NaCl concentration $c*$. We can refer to this condition as the equilibrium response, such that Eq. (2.1) may be specialized to

*π*in cartilage, this study borrowed a concept from the pioneering study of Eisenberg and Grodzinsky [23]. Their approach recognized that, under sufficiently hypertonic conditions, the Donnan osmotic pressure reduces to zero since mobile counterion charges can overwhelmingly shield the fixed charge density

*c*in the tissue ($c*\u226bcF$).

^{F}^{2}By defining a zero-strain reference in this state, the stress required to attain this configuration in a more hypotonic environment can only represent a pressure term. In this study, their concept of a fixed hypertonic reference has been adopted to allow direct measurement of osmotic pressure (illustrated in Fig. 1). The zero-strain reference configuration $\epsilon =0$ of the porous solid matrix was defined as the equilibrium position of the cartilage plug after application of a small tare stress

*σ*under hypertonic conditions; defining that state as the reference configuration produced $\sigma eq(0,2\u2009M)=0$. After stress-relaxation to an equilibrium state, the tissue sample was osmotically loaded by reducing the bath concentration to various desired values $c*$; due to the confining chamber, swelling only occurred in the axial direction. Then, a compressive stress $\sigma 0(c*)$ was prescribed (over and above the tare stress) to bring the cartilage plug back to its zero-strain reference position (the tare configuration). By definition, since $\sigma eqe$ was zero at $\epsilon =0$, the only stress in the tissue was the osmotic pressure

_{t}Therefore, the stress $\sigma 0(c*)$ required to return the tissue to the reference state provided a measurement of the osmotic pressure in situ. Note that Eq. (2.3) requires $\pi (0,2\u2009M)=0$, as the analysis described in this section must be relative to a known value. The error introduced by this assumption was small (Ref. [15], see Sec. 4 for further details).

#### 2.1.2 Osmotic Pressure.

*π*inside articular cartilage may be written as

*R*is the universal gas constant,

*θ*is the absolute temperature, $\Phi $ is the osmotic coefficient inside the tissue, $\Phi *$ is the osmotic coefficient in the external bath, $\gamma \xb1$ is the mean activity coefficient inside the tissue, and $\gamma \xb1*$ is the mean activity coefficient in the external bath. The FCD is a function of strain through the mass balance relation [34]; in a confined compression analysis, it is given by

where $\phi 0w$ is the fluid volume fraction in the reference configuration and $c0F$ is the FCD in the reference configuration; consequently $cF=c0F$ when $\epsilon =0$. In the most general physical chemistry framework, osmotic and activity coefficients (which denote deviations from ideal physicochemical conditions) may be functions of the external bath concentration $c*$, the tissue ion concentrations *c*^{+} and *c*^{–}, and the fixed charge density *c ^{F}*. Of these parameters, values of $\Phi *$ and $\gamma \xb1*$ may be obtained from electrolyte reference tables [38] for a known $c*$; biochemical assays (Sec. 2.5) may be used to determine

*c*. As a result, two unknowns still remain, $\Phi $ and $\gamma \xb1$, though the experiment described in Sec. 2.1.1 is only sufficient to characterize a single unknown.

^{F}*c*in Eq. (2.4) was replaced by an effective fixed charge density $\xi cF$, where $\xi \u22641$ is a scaling parameter. Applying these assumptions to Eq. (2.4) produces

^{F}*c*in Eq. (2.5)) and ion concentrations. Due to the standard assumption of electroneutrality inside the tissue [10,41], and our use of an effective FCD, the concentrations satisfied $c+=\xi cF+c\u2212$ such that an explicit dependence on only one of the two tissue ion concentrations was sufficient. Since cartilage proteoglycans are negatively charged under physiologic pH, we chose the anion concentration

^{F}*c*

^{–}. In particular, under equilibrium conditions,

*c*

^{–}is evaluated as

*X*is defined as

*X*represents the relative effect of the fixed charges. In general,

*c*

^{–}is evaluated by solving the governing equations of triphasic theory; under equilibrium conditions, it is given by Eq. (2.7). Whenever $c*\u226b\xi cF$ (hypertonic conditions), $X\u21920$ and the fixed charges are essentially shielded by the excess mobile ions. Conversely, as $c*\u21920$ (hypotonic conditions), $X\u2192\u221e$ and the tissue behavior is dominated by the presence of the fixed charges. The osmotic coefficient was then defined in analogy to Manning's formula [42] as

Equations (2.6)–(2.9) fully characterize a model for nonideal Donnan osmotic swelling with a single unknown parameter *ξ*, based on the assumption of an effective fixed charge density. *ξ* represents the single free parameter in this model and was determined by fitting to experimental data.

It should be noted that our use of the Manning parameter *ξ* in Eq. (2.9) was constitutive and phenomenological in nature and hence did not imply the microstructural modeling assumptions outlined in the work by Manning [42], especially since that study did not scale the FCD with *ξ* as done here in Eqs. (2.6)–(2.8). Rather, our purpose was to provide a constitutive model for the Donnan osmotic pressure which could be successfully fitted to experimental data on articular cartilage and easily incorporated into theoretical or computational studies that seek to accurately capture osmotic behavior. A consequence of this approach, which is more fully discussed in Sec. 4, is that some physical interpretation (in the sense of physical chemistry) has been lost. For example, there is no experimental evidence showing that $\Phi =\Phi *$ holds in general; however, using experimentally characterized values of $\Phi *$ for various salt solutions would require customized formulations for $\Phi $ that monotonically reduce *π* toward zero in Eq. (2.4) as $c*$ increases. Such solution-specific formulations are neither trivial nor practical for purposes of modeling Donnan osmotic pressure in biomechanics. As the following sections demonstrate, the simplicity and accuracy of our model overrode such concerns.

Finally, it should be noted that *π* is not exactly equal to zero at $c*=2\u2009M$. Therefore, experimental measurements actually represent osmotic pressure differences relative to the unknown osmotic pressure at $c*=2\u2009M$. As discussed further in Sec. 4, a post hoc analysis reveals the error introduced by setting $\pi (0,2\u2009M)\u22480$ is minimal.

### 2.2 Experimental Design.

Three groups were utilized in this study to probe the osmotic pressure over a wide range of FCDs and compare bovine to human tissue. Bovine tissue was split into groups G1 and G2, and human tissue in group H1. Osmotic pressure measurements and confined compression stress-relaxation experiments were performed using a custom-built mechanical loading device [43] to apply axial loading to a final engineering strain $\epsilon =\u2212\epsilon a$. Groups G1 and H1 were loaded to $\epsilon a=0.2$ and group G2 to $\epsilon a=0.1$. Each sample from every group was tested under six different applied bathing concentrations in order of decreasing concentration (see below).

### 2.3 Sample Harvesting and Preparation.

Cylindrical cartilage explants for groups G1 and G2 were harvested from the femoral condyle of knee joints of immature bovine calves (2–3 months old, mixed sex) using a trephine $(\u2205\u20095\u2009mm)$. Explants were placed articular surface down on a sledge microtome (Leica Instruments #SM2400, Nusslock, Germany) equipped with a freezing stage (Physitemp Instruments #BFS-30TC, Clifton, NJ) and embedded and frozen with water-soluble sectioning gel (Thermo-Fisher Scientific #1310APD, Rockford, IL). Any subchondral bone and chondroepiphysis was removed; the sample was then flipped and the superficial and middle zones were removed, producing a final thickness *h*. All samples were stored in 2 M NaCl with a protease inhibitor and phosphate buffers (see below) and frozen at $\u221220\u2009\xb0$C until the day of testing. Resulting specimens were cylindrical plugs of deep zone articular cartilage (G1: $n=10,h=0.70\xb10.18$ mm; G2: $n=8,h=0.76\xb10.02$ mm).

Visually undamaged human cartilage explants (modified Outerbridge grade 0/1; [44,45]) for group H1 were harvested from the patellar groove of knee joints of four fresh-frozen human cadavers obtained from a tissue bank (Anatomy Gifts Registry, Hanover, MD^{3}; two male, two female; aged 65.7±5.4 years) using a trephine $(\u2205\u20095\u2009mm)$. Explants were microtomed in the same manner as bovine tissue. Resulting human specimens were cylindrical plugs primarily containing middle zone cartilage (H1: $n=6,h=0.51\xb10.09$ mm). In this study, multiple samples from the same joint may be considered independent samples, since swelling behavior of plugs were modeled relative to their individual biochemical compositions. All samples were stored in 2 M NaCl with a protease inhibitor and phosphate buffers (see below) and frozen at $\u221220\u2009\xb0$C until the day of testing.

Six NaCl bathing concentrations were prepared for the study, with concentrations $c*$ covering equal increments in log space between hypertonic 2 M and hypotonic 0.001 M (nominally $c*=2,0.6,0.15,0.03,0.006,and\u20090.001$ M NaCl; these concentrations are, respectively, denoted as $M1\u2212M6$ in what follows). The hypertonic solution *M*_{1} contained the prescribed concentration of NaCl supplemented with 0.5 mM ethylenediamine tetra-acetic acid, 35 mM sodium phosphate monobasic, and 0.04% isothialozone-based biocide (Proclin 950, Sigma-Aldrich #46878-U, St. Louis, MO). The full solution information is provided in Table 1. Solutions $M2\u2212M6$ were created based on successive dilutions of *M*_{1}. Use of successive dilutions was necessary to also dilute the buffers and protease inhibitor, as osmotic pressure is a colligative property. Solutions were tested to ensure that pH remained between 6.8 and 7.2. Following each experiment, aliquots of the testing solution were drawn and their osmolarity was tested on a freezing-point osmometer (Model 3250, Advanced Instruments, Norwood, MA). For the remainder of this presentation, $c*$ values reported as concentrations of NaCl are understood to refer to the NaCl/protease inhibitor/buffer solutions described here.

### 2.4 Mechanical Testing.

The experimental protocol described below is summarized in Fig. 2. On the day of testing, samples were thawed in nominal 2 M NaCl $(M1)$ and their thickness and wet weight was measured at this hypertonic configuration. Each sample was then placed inside a stainless steel smooth cylindrical confining chamber $(\u22054.87\u2009mm)$ with impervious sidewalls and bottom surface, machined out of stainless steel. The entire chamber was immersed in a bath containing solution *M*_{1}. A custom mechanical loading device was used to apply axial loading. The top surface of the tissue was loaded in compression with a rigid sintered porous steel indenter $(\u22054.74\u2009mm)$ in series with a 45 N load cell (Model 31, Honeywell Sensing and Control, Charlotte, NC) and LVDT (Fig. 3). Motion was provided by a stepper motor controlled by custom LabVIEW software (LABVIEW 2010, National Instruments Corporation, Austin, TX). Confined compression stress-relaxation tests were performed on each sample while immersed in six different external baths $M1\u2212M6$, in order from highest $(M1)$ to lowest $(M6)$ concentration. The first test, in 2 M NaCl, was used to determine a baseline response for the tissue. In this initial configuration, a 500 kPa preconditioning stress was applied for 60 s and then removed, to securely seat the sample in the confining chamber. Following a 5-min recovery period, a $\sigma t=50\u2009kPa$ tare stress was applied for 15 min (creep equilibrium was attained in ∼5 min), with the reference position for each sample defined as the position at the end of this time [46–49]. By pooling all samples, tare strains relative to the cut thickness were calculated to be 14.8 ± 4.0%, and 17.7 ± 5.7% relative to the now-compressed thickness; these values are in agreement with prior studies demonstrating that 15-20% strain is required for suitable confinement [47–49]. A ramp-and-hold stress-relaxation to a prescribed engineering strain $\epsilon a$ was then applied; after stress equilibration at $\epsilon a$, the porous indenter was lifted, the bathing solution was replaced with the next NaCl concentration, and the sample was allowed to equilibrate to a new swelling state $\epsilon >0$ while still inside the confining chamber (∼1 h; necessary equilibration time estimated based on sample thickness, NaCl diffusivity, and representative cartilage permeability values).

To replace the bathing solution, the confining chamber holding the sample was removed from the bath; the bathing container was then emptied, rinsed in distilled water, dried, and refilled with the next NaCl concentration. To avoid density gradients or poor mixing, any of the previous NaCl solution inside the confining chamber adjacent to the cartilage surface was removed with a tissue paper (Kimwipe^{4}) prior to placing the confining chamber and sample into the new bath. Keeping the sample inside the confining chamber for the duration of the entire experiment was necessary to maintain tight confinement. During the swelling step in the hypotonic range, in rare instances some cartilage plugs swelled unevenly against the side wall of the confining chamber such that they lifted off the bottom and became wedged at an angle in the chamber. When this was observed, the porous indenter tip was used to manually press these samples back down until they contacted the chamber bottom; in such cases, samples were allowed to equilibrate for an additional 15 min.

Prior to beginning the next stress-relaxation test, the porous indenter slowly (∼1 *μ*m/s) applied a prescribed displacement to return the sample to the defined zero-strain reference configuration, where it was again allowed to equilibrate (∼10–20 min). This procedure ensured each stress-relaxation test began at the same reference configuration. The stress $\sigma 0(c*)$ required to return each sample to its nominal zero-strain configuration from its swelled state was measured and used to evaluate the swelling pressure as per Eq. (2.3). In each stress-relaxation test, displacement of the sample was prescribed in the form of a linear ramp at a rate of 0.25 *μ*m/s, and then held constant at the prescribed engineering strain $\epsilon a$ until the stress reached an equilibrium value $\sigma 1(\u2212\epsilon a,c*)$ (∼1 h). For the duration of each stress-relaxation experiment, the bath was covered in plastic wrap only pierced by a hole of sufficient diameter to accept the loading platen, to minimize any evaporation-driven changes in $c*$. All testing for each individual sample was conducted on a single day, and lasted ∼14 h on average. Following testing, each sample was placed in a vacuum desiccator for 72 h and weighed to obtain the tissue dry weight.

To address the potential concerns about the lengthy testing duration, a preliminary study with bovine tissue using the same protocol described here also performed another test in *M*_{1} after completion of testing in solutions $M1\u2212M6$, verifying that the baseline results were maintained. In this preliminary study the *M*_{1} aggregate modulus *H _{A}* was treated as a proxy for tissue degradation over the testing duration. The first and second

*M*

_{1}tests produced aggregate moduli of $HA(1)=0.41\xb10.10$ MPa and $HA(2)=0.38\xb10.17$ MPa, respectively (

*n*=

*10). No significant difference between the groups was seen with a one-tailed paired*

*t*-test (

*p*=

*0.23), confirming that the protocol prevented noticeable tissue degradation during the testing procedure.*

### 2.5 Biochemical Assay.

Dried samples were digested with 0.5 mg/mL proteinase K solution (Fisher Scientific) overnight at $56\xb0$ C, and the FCD was measured by the 1,9-dimethylmethylene blue dye-binding spectrophotometric assay [50], which binds to negative charges in the sample, using shark chondroitin sulfate (CS) as a standard. Given that the charge number of CS is –2, the FCD (normalized to the fluid content of cartilage) could be calculated by also using measurements of cartilage wet and dry weights as described previously [36].

### 2.6 Osmotic Model Fitting.

*s*of a species, the value

*ξ*which provided the best fit over all testing concentrations was determined. We then evaluated a single value $\xi =\xi \xaf$ to model bovine or human cartilage, using the weighted average

_{s}The referential fixed charge density $c0F$ was used here since these osmotic pressure measurements were performed at the defined zero-strain reference state; see Eq. (2.5). Predictions of the mean osmotic pressure response were then made using the mean $c0F$ for the species and the calculated $\xi \xaf$.

### 2.7 Finite Element Modeling.

Finite element curve-fitting with triphasic mixture theory was performed on the transient confined compression stress-relaxation data using the free, open source finite element software febio^{5} [51] to determine the influence of $c*$, *c ^{F}*, and $\epsilon a$ on the mechanical properties of human and bovine cartilage. To account for nonideal osmotic swelling, the osmotic coefficient model given by Eq. (2.9) was implemented in febio as “osm-coef-Manning,” which requires only the specification of

*ξ*and the identification of the co-ion (the anion in this case). For these fits, a more granular value of

*ξ*was adopted for each sample

*s*, such that Eq. (2.3) was satisfied exactly by Eqs. (2.6)–(2.9) for each concentration $c*$ and measured stress $\sigma 0(c*)$; this was necessary to obtain the most accurate fits against experimental data. Since we assumed that $\pi (0,2\u2009M)\u22480$, which is equivalent to assuming $X\u22480$, a meaningful granular value of

*ξ*is not available at that concentration according to Eq. (2.8). Therefore, at $2\u2009M$, the fits were performed as though cartilage was essentially a biphasic material (triphasic with zero FCD), in agreement with the assumption $\pi (0,2\u2009M)\u22480$.

*c*

^{–}representing the co-ion. The neo-Hookean strain energy density was given by

*G*and

*λ*are material coefficients, $J=detF$ is the Jacobian of the deformation gradient

**F**, and $I1=tr\u2009C$, where $C=FT\xb7F$ is the right Cauchy–Green tensor [52]. The material coefficients are related to Young's modulus

*E*and Poisson's ratio

*ν*through $E=2G(1+\nu )$ and $\lambda =2G\nu /(1\u22122\nu )$. Holmes-Mow strain-dependent permeability, which considers the reduction in permeability as the interstitial pores in cartilage collapse upon loading, was employed based on preliminary data and experimental results in the literature [46,53,54], where we adopted the form

The two adjustable parameters in Eq. (2.12) are the zero-strain permeability *k*_{0} (same units as permeability *k*) and the exponential parameter *m* (unitless). For curve-fitting, the material properties to be fitted were taken to be *E*, *k*_{0}, and *m*; *ν* = 0 was prescribed due to the use of a one-dimensional confined compression analysis and the assumption that the pore volume can be reduced with negligible lateral expansion of the porous matrix. Consequently, the fitted property *E* also represents the cartilage zero-strain aggregate modulus *H _{A}*.

^{2}/s and $D0\u2212=D\u2212=0.0016$ mm

^{2}/s [55,56]. This assumption corresponds to neglecting hindrance of solute flux by the porous solid matrix; however, the effective solute diffusivities inside the tissue could be affected by charge interactions due to the FCD [12,57]. In febio, boundary conditions on the tissue are prescribed on the solid displacement, the effective fluid pressure $p\u0303$ and the effective solute concentration $c\u0303\alpha $ [34]. In this confined compression model, in addition to displacement boundary conditions that enforce the confining conditions, boundary conditions needed to be prescribed at the free-draining porous indenter. Based on the assumptions of Sec. 2.4, the boundary conditions were

^{6}

Due to our assumption that $\Phi *=\Phi $ in this constitutive model for the osmotic pressure, and since $\Phi $ varies over space and time according to Eqs. (2.8) and (2.9), a specialized boundary condition, called “matching_osm_coef,” was implemented in febio to enforce this constraint.

The finite element mesh employed in the curve-fitting procedure used 30 isoparametric hexahedral elements in the axial direction, biased to be thinner toward the free-draining indenter. Only the axial direction needed to use mesh refinement, since confined compression is a one-dimensional analysis. Curve-fitting was performed through febio's Parameter Optimization module, which uses a Levenberg–Marquardt algorithm. To avoid local minima and nonuniqueness of converged solutions, ten sets of randomized initial guesses were used as input values to the algorithm for each fit, and the converged material parameters producing the lowest error were selected. For each concentration $c*$, ten initial guesses were generated for each material parameter $(HA,\u2009k0,\u2009m)$ based on prior studies in the literature [15,35,46,47,54] and the spread of results from preliminary studies. The average runtime for a model with specified material parameters was 90–180 s; full parameter optimization required 15–90 min for each set of initial guesses, depending on the number of model runs required.

### 2.8 Statistical Analyses.

Shapiro–Wilk tests with a significance level of *p *=* *0.05 were performed on all measured and fitted values, to determine whether parametric or nonparametric statistics were appropriate. If the Shapiro–Wilk test rejected the null hypothesis of a normal distribution, the command FindDistribution was used in Mathematica 11 (Wolfram Research, Champaign, IL) to determine the statistical distribution which best represented the data. If a lognormal distribution was indicated, the Shapiro–Wilk test was then performed on the log-transformed data.

Experiments were performed in six different concentrations $c*$ for groups G1, G2, and H1; additionally, groups G1 and G2 were tested at different applied strains. For normally distributed data in groups G1 and G2, two-way analysis of variance (ANOVA) ($\alpha =0.05$) was performed to seek significant main effects of concentration (between $M1\u2212M6$ tests) and applied strain ($\epsilon a=0.2$ and 0.1) and all possible interactions, with Tukey's HSD post hoc testing of the means implemented when significance was detected (*p *<* *0.05). If no main effect of applied strain was seen, data from G1 and G2 were pooled for the given analysis. For data which did not follow a normal or lognormal distribution, two-way ANOVA was still performed, due to the robustness of ANOVA to non-normality [58,59]; however, these results should be interpreted with more care. As human data was only obtained at a single strain value, one-way ANOVA for the factor $c*$ was employed to analyze group H1, with the same caveat regarding normality.

In the presentation of results, lower-case letters convey main effects of concentration, according to the appropriate ANOVA. Letters shared between data points indicate those points which are not significantly different (*p *>* *0.05). ANOVA $p\u2212$values for main effects of concentration $c*$ and applied strain $\epsilon a$ are directly reported in each plot. Interactions $\epsilon a\xd7c*$ are reported in figure captions, where applicable. ANOVA results only communicate results within a single plot, hence shared letters between plots have no meaning. All presented numerical values represent mean ± one standard deviation.

## 3 Results

### 3.1 Experimental Measurements.

Representative sets of experimental stress-relaxation curves are shown in Fig. 4 (symbols) for both bovine (G1) and human (H1) samples. Every curve in Fig. 4 is produced by the same sample under the same loading profile. Differences are only due to a change in external bathing concentration $c*$. For each sample, $\sigma 0(c*)$ was directly measured from the first data point (*σ* at time *t *=* *0) in each individual stress-relaxation curve; the curves start at different points due to the osmotic pressure, with the initial stress equal to zero at $c*=2$ M. The equilibrium stresses $\sigma 1(c*)$ were similarly obtained from the final points of each stress-relaxation curve. Though the magnitudes differed between G1, G2, and H1, the set of curves for every sample were substantially similar qualitatively. From hypertonic to approximately isotonic conditions, there is an inflection point in the stress during the ramp phase, similar to that predicted by Holmes et al. [46] (this is more easily seen for $M1\u2212M3$ in Fig. 4(a)). As $c*$ becomes progressively more hypotonic the stress response becomes concave downward. By normalizing all stress-versus-time relaxation curves for a single sample to their respective equilibrium stress $\sigma 1(c*)$ and ramp time, it becomes evident that the relaxation time increases with increasing $c*$ (Fig. 5). Normalized curves between samples were strikingly similar. A related finding was that the ratio of normalized peak to equilibrium stresses increased monotonically from a plateau value as $c*$ increased (Fig. 5), with very little scatter in the data. Two-way ANOVA revealed statistically significant effects of concentration and strain for groups G1 and G2; a corresponding one-way ANOVA indicated an effect of concentration for group H1.

The osmolarity of each solution $M1\u2212M6$ was directly measured using a freezing point osmometer. Measured values are reported in Table 2, showing deviations from nominal $c*$ values. The deviations in $M2\u2212M6$ are likely associated with the successive dilution process. Deviation in *M*_{1} from $c*=2$ M may have originated from the use of insufficiently precise measuring equipment, as large quantities of this solution were prepared at a given time.

M_{1} (2 M) | M_{2} (0.6 M) | M_{3} (0.15 M) | M_{4} (0.03 M) | M_{5} (0.006 M) | M_{6} (0.001 M) |
---|---|---|---|---|---|

1.836 ± 0.013 | 0.543 ± 0.005 | 0.138 ± 0.0005 | 0.027 ± 0.0003 | 0.004 ± 0.0003 | 0.001 ± 0.0 |

M_{1} (2 M) | M_{2} (0.6 M) | M_{3} (0.15 M) | M_{4} (0.03 M) | M_{5} (0.006 M) | M_{6} (0.001 M) |
---|---|---|---|---|---|

1.836 ± 0.013 | 0.543 ± 0.005 | 0.138 ± 0.0005 | 0.027 ± 0.0003 | 0.004 ± 0.0003 | 0.001 ± 0.0 |

Values are reported in units of M; nominal values for each measurement are displayed in the top row. The measured values were calculated from the osmolarity by assuming only NaCl was present. No standard deviation was associated with *M*_{6}, as the measurements were at the limit of the device resolution (1 mOsm).

Measured osmotic pressures satisfied the hypothesis of a normal distribution, according to Shapiro–Wilk testing. No significant differences in measured osmotic pressure between G1 and G2 were detected by two-way ANOVA $(\epsilon a:\u2009p=0.26)$, thus the bovine data were pooled. The osmotic pressure for the pooled bovine data (Fig. 6(a)) monotonically increased as $c*$ decreased, reaching a steady value when $c*\u22640.006$ M. The osmotic pressure in H1 (Fig. 6(b)) also rose monotonically, yet the trend was not as smooth, likely owing to greater variability in human tissue samples and the smaller sample size. Nevertheless, the increase in osmotic pressure in H1 was statistically significant. The measured osmotic pressure in H1 only achieved about two-thirds of the values reached by the pooled G1–G2 samples.

The measured referential fixed charge density in G1 was $c0F=344\xb134\u2009mEq/L$; according to a two-tailed homoscedastic t-test, this was not significantly different than $c0F=315\xb134\u2009mEq/L$ measured for G2 ($p=0.12)$, further supporting the decision to combine both groups in Fig. 6(a). Human samples in H1 displayed a lower FCD with greater variability ($c0F=231\xb159\u2009mEq/L$), a possible explanation for the greater variance seen in the H1 osmotic pressure data in Fig. 6(b).

Referential fluid volume fractions $\phi 0w$ for groups G1, G2, and H1 were calculated from wet and dry weights and determined to be $\phi 0w=0.743\xb10.009,\u2009\phi 0w=0.733\xb10.010$, and $\phi 0w=0.698\xb10.015$, respectively. Following compression, Eq. (2.5) gives $cF=471\xb149\u2009mEq/L$ for G1 compared to $cF=365\xb141\u2009mEq/L$ in G2, a significant difference $(p<0.001)$. In H1, Eq. (2.5) reports $cF=324\xb184\u2009mEq/L$.

The equilibrium stress $\sigma 1(c*)$ at the end of stress-relaxation (Fig. 7) followed the same trend as the zero-strain osmotic pressure data. Equilibrium stress in G1 was larger than G2, which was expected due to differences in applied strain $\epsilon a$ between the groups. Equilibrium stresses in human tissues (Fig. 7(b)) were comparable to those observed in G1 with the same $\epsilon a$ (Fig. 7(a)), despite differences up to 100 kPa in osmotic pressure at zero strain. All three groups again displayed similar trends, showing a log-sigmoidal increase in equilibrium stress as $c*$ decreased, plateauing when $c*\u22640.006$ M.

### 3.2 Osmotic Pressure.

By adjusting the parameter *ξ*, the nonideal Donnan model described by Eqs. (2.6)–(2.9) was able to provide an excellent fit to osmotic swelling data from both bovine and human cartilage. Sample-specific values *ξ _{s}* ranged from 0.27 to 0.77. The first modeling efforts considered pooled bovine (G1–G2) and human (H1) groups separately, with the species-averaged parameter $\xi \xaf$ calculated according to Eq. (2.10). This parameter was found to be $\xi \xaf=0.472$ for bovine cartilage and $\xi \xaf=0.412$ for human cartilage. Using $\xi \xaf$ and mean $c0F$ values to model the osmotic pressure of G1–G2 and H1 produced good agreement ($R2=0.89$ for bovine and $R2=0.76$ for human; Fig. 8, dashed curves).

Combining all groups (G1, G2, H1) produced $\xi \xaf=0.461$, and applying this $\xi \xaf$ to both bovine and human cartilage, together with the appropriate average $c0F$, also produced good predictions of the osmotic pressures ($R2=0.89$ for bovine and $R2=0.73$ for human; Fig. 8, solid curves). These results suggest the simple framework of Eqs. (2.6)–(2.9), together with the single parameter $\xi \xaf=0.461$, may be used to accurately simulate osmotic swelling of articular cartilage for both species.

Granular values of *ξ* ranged from 0.107 to 1.00. Sample-specific values of the osmotic coefficient $\Phi $ from all tests were extracted by using these granular values in Eq. (2.9) (symbols in Fig. 9). These values were in close agreement with the smooth response provided by Eq. (2.9) when using $\xi \xaf=0.461$ ($R2=0.64$; solid curve in Fig. 9). That curve is bounded by the values $\Phi |X\u21920=1$ and $\Phi |X\u2192\u221e=1\u221212\xi \xaf$.

### 3.3 Material Properties.

Material properties *H _{A}*,

*k*

_{0}, and

*m*were extracted by curve-fitting stress-relaxation data sets, such as those shown in Fig. 4, using the granular values of

*ξ*for each bath concentration, and the corresponding effective FCD, $\xi c0F$. Curve-fits for representative bovine (G1) and human (H1) samples are shown in Fig. 4. Similarly good fits were achieved for tests on all samples. A total of 1440 curve-fits were performed in this study (24 total samples × 6 tests per sample × 10 sets of initial guesses ${HA,k0,m}$ to avoid local minima). Average

*R*

^{2}values were 0.986±0.007 for G1, 0.975±0.013 for G2, and 0.977±0.015 for H1. Tables 3–5 show the goodness of fit values broken down by concentration.

G1 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.436 ± 0.073 | 1.272 ± 0.310 | 8.58 ± 0.49 | 0.987 ± 0.005 |

0.54 | 0.542 ± 0.112 | 1.597 ± 0.457 | 8.23 ± 1.05 | 0.985 ± 0.009 |

0.138 | 0.785 ± 0.192 | 1.232 ± 0.337 | 8.34 ± 1.02 | 0.987 ± 0.008 |

0.027 | 1.080 ± 0.211 | 0.842 ± 0.243 | 7.89 ± 0.89 | 0.981 ± 0.009 |

0.004 | 1.093 ± 0.207 | 0.633 ± 0.150 | 6.66 ± 0.73 | 0.986 ± 0.004 |

0.001 | 1.105 ± 0.176 | 0.580 ± 0.095 | 6.17 ± 0.46 | 0.988 ± 0.003 |

G1 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.436 ± 0.073 | 1.272 ± 0.310 | 8.58 ± 0.49 | 0.987 ± 0.005 |

0.54 | 0.542 ± 0.112 | 1.597 ± 0.457 | 8.23 ± 1.05 | 0.985 ± 0.009 |

0.138 | 0.785 ± 0.192 | 1.232 ± 0.337 | 8.34 ± 1.02 | 0.987 ± 0.008 |

0.027 | 1.080 ± 0.211 | 0.842 ± 0.243 | 7.89 ± 0.89 | 0.981 ± 0.009 |

0.004 | 1.093 ± 0.207 | 0.633 ± 0.150 | 6.66 ± 0.73 | 0.986 ± 0.004 |

0.001 | 1.105 ± 0.176 | 0.580 ± 0.095 | 6.17 ± 0.46 | 0.988 ± 0.003 |

All values reported are mean ± one standard deviation. To obtain the true values of *k*_{0}, multiply the table values by $10\u22123$.

G2 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.400 ± 0.140 | 1.332 ± 0.339 | 8.10 ± 1.37 | 0.982 ± 0.009 |

0.54 | 0.518 ± 0.130 | 1.874 ± 0.616 | 9.34 ± 2.12 | 0.979 ± 0.011 |

0.138 | 0.884 ± 0.372 | 1.789 ± 0.745 | 10.43 ± 3.04 | 0.971 ± 0.015 |

0.027 | 1.258 ± 0.434 | 1.619 ± 0.737 | 11.52 ± 2.71 | 0.973 ± 0.013 |

0.004 | 1.386 ± 0.403 | 1.144 ± 0.476 | 10.08 ± 2.85 | 0.967 ± 0.012 |

0.001 | 1.496 ± 0.576 | 1.102 ± 0.392 | 9.92 ± 2.98 | 0.979 ± 0.009 |

G2 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.400 ± 0.140 | 1.332 ± 0.339 | 8.10 ± 1.37 | 0.982 ± 0.009 |

0.54 | 0.518 ± 0.130 | 1.874 ± 0.616 | 9.34 ± 2.12 | 0.979 ± 0.011 |

0.138 | 0.884 ± 0.372 | 1.789 ± 0.745 | 10.43 ± 3.04 | 0.971 ± 0.015 |

0.027 | 1.258 ± 0.434 | 1.619 ± 0.737 | 11.52 ± 2.71 | 0.973 ± 0.013 |

0.004 | 1.386 ± 0.403 | 1.144 ± 0.476 | 10.08 ± 2.85 | 0.967 ± 0.012 |

0.001 | 1.496 ± 0.576 | 1.102 ± 0.392 | 9.92 ± 2.98 | 0.979 ± 0.009 |

All values reported are mean ± one standard deviation. To obtain the true values of *k*_{0}, multiply the table values by $10\u22123$.

H1 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.499±0.208 | 0.849±0.276 | 8.98±1.80 | 0.969±0.021 |

0.54 | 0.575±0.210 | 1.101±0.212 | 8.79±0.74 | 0.976±0.014 |

0.138 | 0.885±0.321 | 1.318±0.673 | 9.88±2.50 | 0.976±0.009 |

0.027 | 1.169±0.500 | 0.753±0.249 | 7.92±2.46 | 0.974±0.018 |

0.004 | 1.429±0.496 | 0.689±0.304 | 7.63±2.86 | 0.983±0.007 |

0.001 | 1.597±0.455 | 0.863±0.379 | 8.86±2.98 | 0.985±0.006 |

H1 | ||||
---|---|---|---|---|

$c*$ (M) | H (MPa)_{A} | $k0\xd7103$ (mm^{4}/N$\xb7$s) | m | R^{2} |

1.84 | 0.499±0.208 | 0.849±0.276 | 8.98±1.80 | 0.969±0.021 |

0.54 | 0.575±0.210 | 1.101±0.212 | 8.79±0.74 | 0.976±0.014 |

0.138 | 0.885±0.321 | 1.318±0.673 | 9.88±2.50 | 0.976±0.009 |

0.027 | 1.169±0.500 | 0.753±0.249 | 7.92±2.46 | 0.974±0.018 |

0.004 | 1.429±0.496 | 0.689±0.304 | 7.63±2.86 | 0.983±0.007 |

0.001 | 1.597±0.455 | 0.863±0.379 | 8.86±2.98 | 0.985±0.006 |

All values reported are mean±one standard deviation. To obtain the true values of *k*_{0}, multiply the table values by $10\u22123$.

*H*,

_{A}*k*

_{0}, and

*m*for all groups, along with

*R*

^{2}values from febio curve-fitting. The aggregate modulus

*H*was lognormally distributed for groups G1–G2 (Shapiro-Wilk on $lnHA$:

_{A}*p*=

*0.21) and normally distributed for H1 (Shapiro–Wilk:*

*p*=

*0.12), thus allowing the use of parametric statistics. Two-way ANOVA detected no statistical significance for*

*H*between groups G1 and G2 $(\epsilon a:\u2009p=0.26)$; accordingly, bovine

_{A}*H*data were pooled. Similar to trends seen for osmotic pressure, equilibrium stress, and osmotic coefficient, the

_{A}*H*data were strongly concentration-dependent for both bovine $(p<0.0001)$ and human $(p<0.001)$ groups, with a log-sigmoidal shape (Fig. 10). A model of the form

_{A}was fit to each *H _{A}* dataset (pooled G1-G2 and H1) to extract a continuous function $HA(X)$ with units of MPa. In Eq. (3.1), the modulus as $X\u2192\u221e$ is denoted $HA,\u221e$, whereas $HA,\delta $ represents the increase in modulus relative to $HA|X\u21920$. To maintain consistency with the proposed model for nonideal Donnan osmotic swelling,

*X*was calculated at $\epsilon =0$ using $\xi \xaf=0.461$ in the model of Eq. (3.1). Curve-fitting parameters for pooled bovine (G1-G2) and human (H1) groups are reported in Table 6. Equally good fits were achieved for pooled bovine and human data. The fitted model displayed a qualitatively similar shape between groups, although the plateau as $X\u2192\u221e$ was not as stark in H1. This lack of a plateau was in large part due to the exclusion of a single outlier in one human sample at

*M*

_{6}; see Section 4 for more detail.

Group | a | b | $HA,\delta $ (MPa) | $HA,\u221e$ (MPa) | R^{2} |
---|---|---|---|---|---|

G1–G2 | 1.305 | 0.738 | 0.940 | 1.246 | 0.914 |

H1 | 1.902 | 0.523 | 1.160 | 1.553 | 0.896 |

Group | a | b | $HA,\delta $ (MPa) | $HA,\u221e$ (MPa) | R^{2} |
---|---|---|---|---|---|

G1–G2 | 1.305 | 0.738 | 0.940 | 1.246 | 0.914 |

H1 | 1.902 | 0.523 | 1.160 | 1.553 | 0.896 |

Parameters *a* and *b* are unitless.

Permeability parameters *k*_{0} and *m* displayed concentration-dependence in bovine groups G1 and G2 ($k0:\u2009p<0.0001$; $m:\u2009p=0.02$); when comparing G1 and G2, both *k*_{0} and *m* revealed strain-dependence $(k0:\u2009p<0.0001;\u2009m:\u2009p<0.0001)$ (Figs. 11(a) and 11(b)). The parameter *k*_{0} obeyed neither a normal (Shapiro–Wilk: $p<0.0001)$ nor a lognormal (*p *=* *0.006) distribution for bovine (G1-G2) groups. However, a lognormal distribution provided the best fit. In contrast, *k*_{0} followed a lognormal distribution for the human (H1) group (Shapiro–Wilk on $lnk0$: *p *=* *0.24). To analyze *k*_{0}, two-way ANOVA was performed on the log-transformed values for both bovine (G1–G2) and human (H1) groups, though the bovine results should be interpreted in light of the Shapiro–Wilk testing outcome. Similarly, the parameter *m* did not obey either a normal distribution for G1–G2 (Shapiro–Wilk: *p *<* *0.0001), nor was it lognormally distributed $(p=0.041)$. However, a lognormal distribution provided an adequate fit to the data, thus two-way ANOVA was still performed on log-transformed *m* data, though results should be interpreted accordingly. Of note is that a significant interaction effect between strain and concentration was found for *m* in G1-G2 ($\epsilon a\xd7c*:\u2009p=0.004$). For group H1, *m* was normally distributed (Shapiro–Wilk: *p *=* *0.66). In general, both permeability parameters for G1 and G2 were relatively lower at hypertonic and hypotonic concentrations, rising at isotonic concentrations. Human data showed similar trends (Figs. 11(c) and 11(d)), but no statistical significance was detected across concentrations $c*$ using one-way ANOVA $(k0:\u2009p=0.10;\u2009m:\u2009p=0.70)$. The magnitudes of *k*_{0} and *m* for H1 closely matched those of G1-G2. Species-averaged permeability values were $k0=(1.222\xb10.598)\xd710\u22123\u2009mm4/N\xb7s$ and $m=8.61\xb12.32$ for bovine cartilage, and $k0=(0.931\xb10.440)\xd710\u22123\u2009mm4/N\xb7s$ and $m=8.67\xb12.45$ for human cartilage.

## 4 Discussion

This study presents direct measurements of the osmotic swelling pressure in bovine and human articular cartilage explants as a function of bath concentration $c*$, as well as constitutive models for material parameters within the context of triphasic theory [10], demonstrating the dependence of the mechanical response of cartilage on the parameter *X* in Eq. (2.8). The Donnan osmotic pressure was modeled assuming a specialized form of nonideal behavior, where the measured fixed charge density was effectively attenuated by a scaling factor *ξ*, Eq. (2.6). Borrowing from a classical formulation in the physical chemistry literature [42], we adopted the functional form of Eq. (2.9) to model the osmotic coefficient $\Phi $ as a function of *ξ*, the fixed charge density *c ^{F}*, and the anion concentration

*c*

^{–}, accounting for tissue strain via Eq. (2.5). A major outcome of these modeling efforts was the demonstration that a single average parameter $\xi \xaf=0.461$ sufficed to account for the experimental variation in osmotic pressure with FCD, ion concentrations, and tissue strain, for both immature bovine and mature human tissue (Fig. 8). Though minor differences existed between bovine and human cartilage, as revealed by species-specific fitting (Fig. 8, dashed lines), the parameter $\xi \xaf=0.461$ nevertheless was able to capture the experimental results to a similarly high degree of accuracy (Fig. 8, solid lines), particularly given the greater variability in human cartilage. This is beneficial from a modeling perspective, as species-, joint-, or zone-specific values may not always be available. Furthermore, results demonstrated that the material properties of cartilage are a function of

*X*, even within a triphasic framework which explicitly accounts for osmotic contributions (Figs. 10 and 11). The emergence of additional concentration-dependent stiffening not captured by the Donnan model suggests the presence of other coexisting charge effects, which we suggest include Poisson–Boltzmann electrostatic interactions. In addition, cartilage displayed strongly nonideal swelling in the hypotonic range ($\Phi \u21920.764$ as $X\u2192\u221e$ in Fig. 9), in contrast to the common assumption of ideality in the majority of triphasic modeling studies [24,34,57].

The novel constitutive model for predicting the osmotic pressure in cartilage, presented in Sec. 2.1.2, was able to successfully fit the experimentally measured osmotic pressure for all samples in this study, using a single sample-specific parameter *ξ _{s}* for all tested concentrations. The resulting observed dependence of the osmotic coefficient $\Phi $ on concentration (Fig. 9) was qualitatively similar to results from the Poisson–Boltzmann cell model for cartilage osmotic pressure presented by Buschmann and Grodzinsky [15], indicating that these distinct modeling approaches are able to capture the same underlying physical phenomenon. A key finding of this study was that the mean osmotic pressure response for both bovine and human cartilage could be predicted to a high degree of accuracy using a single parameter $\xi \xaf=0.461$ (Fig. 8), calculated from a weighted average of sample-specific

*ξ*values for both species. Physically, this result suggests that osmotic swelling in human and bovine cartilage obeys substantially identical behavior, governed by the ratio of fixed charges to co-ions. From a modeling perspective, when no direct measurements of osmotic pressure are available for a particular cartilage sample or species, Fig. 8 supports the use of $\xi \xaf=0.461$ and the constitutive model of Sec. 2.1.2 to simulate swelling behavior in computational models of articular cartilage. To facilitate modeling endeavors, this constitutive relation was implemented in the free, open-source febio software (starting with febio 3.0).

_{s}A closer examination of Fig. 8 reveals the sensitivity of the osmotic pressure response to $\xi \xaf$ over the range of external bath concentrations, since curves are displayed for two different values of this parameter for both bovine and human results. The sensitivity of the response to *ξ* is further evident if one considers that specimen-specific granular values at each concentration were fitted such that the predicted osmotic pressure matched the corresponding experimental result exactly. Since Fig. 8 displays mean and standard deviations of experimental measurements at each concentration, one can infer that the granular range $0.107\u22121.00$ was needed to match that data scatter exactly. It should be noted that the model in this study was applied to cartilage with FCDs ranging from $150\u2009to\u2009500\u2009\u2009mEq/L$; though there is no reason to believe the model cannot apply to any FCD, it has not been validated outside this range. In particular, whether a single value of $\xi \xaf$ can reliably predict the osmotic response of cartilage whose GAGs have been digested enzymatically remains to be determined from future experiments. Importantly, there is no theoretical impediment to using this value of $\xi \xaf$ for such samples, since our model equations show that the osmotic coefficient approaches unity and the osmotic pressure approaches zero as $c0F\u21920$, regardless of the value of $\xi \xaf$.

This study also examined the transient stress-relaxation response of these cartilage samples under confined compression, at six distinct bath concentrations $c*$, providing a richer dataset than prior related experimental studies. Analyzing the six stress-relaxation curves for each sample showed characteristic trends that were preserved across all groups. Ramp phases of the stress-relaxation curves showed the transient response associated with strain-dependent permeability at higher $c*$, while this effect vanished as $c*$ became more hypotonic (Fig. 4). This finding was supported by the normalized curves, which showed a sharp decrease and eventual plateau in the peak-to-equilibrium stress ratio as $c*$ decreased from *M*_{1} to *M*_{6} (Fig. 5). In an earlier study, an asymptotic analysis of stress-relaxation had linked this ratio to strain-dependent permeability and the tissue's time constant [46]. Visually, it is apparent from Fig. 5 that the relaxation time increases with $c*$, indicating decreasing interstitial fluid diffusivity with greater salt concentrations.

Curve-fitting the transient response provided material properties of cartilage across a wide range of concentrations using the framework of triphasic theory. These material properties showed significant concentration-dependence (Figs. 10 and 11). This dependence is allowed within the general triphasic framework, since triphasic theory includes strain and solute concentrations as state variables. It demonstrated that the proteoglycans not only contribute to the osmotic pressure via their electric charge (i.e., the Donnan effect, which attracts a higher concentration of ions in the tissue than in the bath), they also provided electrostatic (but nonosmotic) contributions to the aggregate modulus *H _{A}* and, possibly, the zero-strain permeability

*k*

_{0}and exponent

*m*(for bovine cartilage). The statistically strong monotonic increase of

*H*with decreasing $c*$ (equivalently, increasing

_{A}*X*; Fig. 10) is consistent with the expectation that proteoglycans become stiffer as the ion concentrations inside the tissue decrease, thereby reducing the charge-shielding effect and producing greater charge-to-charge repulsion forces between glycosaminoglycan side chains, just as predicted by the Poisson–Boltzmann model [15] and experimental measurements [31,33,60]. The Donnan model used in this study characterized a pressure term arising due to a concentration imbalance between the interstitial fluid in cartilage and an external bath, but did not account for tensorial stresses that may have arisen due to electrostatic repulsion. The strong relationship between aggregate modulus and concentration suggests the coexistence of Donnan osmotic pressure and additional charge effects. In particular, the aggregate modulus of both bovine and human tissue increased almost exactly threefold as $c*$ decreased from hypertonic M1 to hypotonic M6 (bovine: $HA=0.420\xb10.109$ MPa to 1.266 ± 0.438 MPa; human: $HA=0.499\xb10.208$ MPa to 1.597 ± 0.455 MPa). The variation of

*H*with the normalized concentration

_{A}*X*followed a smooth log-sigmoidal trend in both species, captured with the simple expression in Eq. (3.1). The classical study by Eisenberg and Grodzinsky [23] previously reported this effect in bovine cartilage, and our results are remarkably similar to theirs.

The variations of *k*_{0} and *m* with $c*$ in bovine samples (Figs. 11(a) and 11(b)) were not monotonic. Furthermore, both parameters significantly varied with $\epsilon a$. In contrast to *H _{A}*, there is no compelling physical argument that would explain the variation of

*k*

_{0}and

*m*with $c*$ or $\epsilon a$, therefore we believe that this dependence may have been an artifact of simplifying constitutive assumptions adopted in our triphasic analysis of the transient response of cartilage, namely, the assumption that ion diffusivities $D0\alpha $ and $D\alpha $ and solubility $\kappa \alpha $ [34] were constant. Experimental evidence exists to suggest that these properties may be dependent on ion concentrations and tissue strain, though they have not been fully characterized for cartilage [55,56,61,62]. These simplifying assumptions directly influence the effective hydraulic permeability $k\u0303$ of the tissue, as described in equations presented in our earlier study [34]. Consequently, while the permeability parameters reported in this study were suitable for successfully fitting the transient stress-relaxation response, they may not have reflected an “intrinsic” or physically meaningful set of values. Therefore, for computational studies of cartilage mechanics in a triphasic framework, we recommend using the species-averaged values, $k0=(1.222\xb10.598)\xd710\u22123\u2009mm4/N\xb7s$ and $m=8.61\xb12.32$ for bovine cartilage, and $k0=(0.931\xb10.440)\xd710\u22123\u2009mm4/N\xb7s$ and $m=8.67\xb12.45$ for human cartilage.

The direct measurements of the osmotic pressure within articular cartilage reported in this study may be compared to measurements of the osmotic pressure of isolated GAG solutions, often used in the literature to indirectly estimate osmotic pressure in cartilage. In an earlier study from our group [30], the osmotic pressure of CS solutions mixed with NaCl was measured for various concentrations of CS and $c*$. The results of that study demonstrated that GAG solutions exhibit a very significant osmotic pressure even at $c*=2\u2009M$; this response was attributed to the phenomenon of configurational entropy [19,29,63]. Therefore, a direct comparison of *π* from this study against Fig. 4 of Chahine et al. [30] shows that the osmotic pressure of GAG solutions considerably overestimates the osmotic pressure measured in cartilage. However, if the osmotic pressure of GAG solutions from configurational entropy (i.e., measurements at $2\u2009M$) is subtracted from the pressures measured at lower concentrations, we find much better agreement between cartilage and GAG solutions, as presented in Fig. 12. These results suggest that that the differential osmotic pressure $\pi (c*)\u2212\pi (2\u2009M)$ measured from GAG solutions provides a reasonable estimate of the osmotic pressure within cartilage; they also suggest that configurational entropy has little influence on the osmotic pressure within cartilage, consistent with the understanding that proteoglycan macromolecules are tightly enmeshed within the collagen matrix, constraining random configurational changes.

There were some experimental limitations to this study. First, the assumption $\pi (0,2\u2009M)\u22480$ adopted in Section 2.1.1 represents an approximation of hypertonic conditions as $c*$ becomes much greater than $\xi cF$ in Eq. (2.6). From a post hoc theoretical analysis, using the overall average value $\xi \xaf=0.461$ in the constitutive model for the osmotic pressure shows that $\pi (0,2\u2009M)$ would have averaged $0.008\xb10.002\u2009MPa$ for bovine samples and $0.004\xb10.002\u2009MPa$ for human samples. Compared to the range of measured osmotic pressures in Fig. 6, we conclude that neglecting this value was reasonable.

Another experimental limitation occurred during testing, when cartilage swelling in rare instances caused the plug to lift off the bottom of the chamber and become wedged at an angle, requiring it to be pressed back down. Of the 144 tests performed in this study (24 samples × 6 tested concentrations), the measured osmotic pressure was significantly larger than predicted from ideal Donnan law in three cases (one hypotonic concentration in each of three different samples) and those measures appeared to be outliers when compared to the smooth trend seen in the remaining five measurements for those samples. In retrospect, we believe that in those three cases the sample was not pressed flush against the bottom of the confining chamber after becoming wedged; those three cases were excluded from the results and analysis. The exclusion of these three data points (one each at concentration *M*_{5} in G2, concentration *M*_{6} in G2, and concentration *M*_{6} in H1) has affected the reported results. In particular, the sample in group H1 which had the *M*_{6} test excluded exhibited the lowest modulus of any H1 samples at each concentration; including this sample's material properties in $M1\u2212M5$ results but removing it from *M*_{6} allows the modulus data for *M*_{6} in group H1 to be slightly skewed upward.

A third potential concern was the lengthy duration of testing, and the repeated mechanical tests performed on each sample. We placed great priority on arresting potential tissue degradation, using an isothiazolone-based biocide along with a protease inhibitor and buffers in the bathing solutions, based on our previous study which demonstrated this combination could maintain devitalized tissue integrity over 28 days at 37 °C [64]. Nevertheless, since it was still possible that damage occurred due to swelling within a confining chamber or repeated loading events, we investigated this possibility in a preliminary study, verifying that the baseline results were maintained.

Another potential limitation concerns cartilage heterogeneity. The depth-dependence of compositional and material properties in both bovine and human cartilage has been well established [36,65–67], with distinct zonal variations seen. However, in our experiments deep zone bovine cartilage and middle zone human cartilage were used, due to the higher FCDs present in these regions; thus all samples in this study were harvested from a single zone and may be considered approximately homogeneous. The model proposed in Sec. 2.1.2 can be incorporated into a depth-dependent study, where zonal variations in *ε* and *c ^{F}* may influence $\Phi $ and

*π*through the tissue depth; the concentration-dependent aggregate modulus

*H*is likely to be a function of depth as well. However, performing such a characterization would be challenging and require experimental validation with digital image correlation [33,67–71], and was beyond the scope of this work.

_{A}Finally, it would have been beneficial to test more than six human samples to strengthen the findings in that species, though external circumstances prevented the completion of those tests.

The use of an effective FCD, given by $\xi cF$ in this study, is consistent with the counterion condensation theory proposed by Manning [42], whereby counterions may condense on the proteoglycans such that their charge density is effectively reduced. However, our modeling approach goes further than Manning's theory, which defined *X* as $cF/c\u2212$ (in contrast to our Eq. (2.8)) and allowed cases where $\xi <1$ and $\xi >1$ (providing an additional formula for $\Phi $ in Eq. (2.9) for that latter case). Manning's theory was only concerned with the osmotic coefficient within a polyelectrolyte solution, whereas the triphasic theory seeks to model osmotic effects relative to an external bath environment. Whereas values of $\Phi *$ and $\gamma \xb1*$ may be found in reference tables for some electrolyte solutions [38] and could be substituted into Wells' extension of Manning's theory [72], we opted for a simpler modeling approach assuming that $\Phi =\Phi *$ and $\gamma \xb1=\gamma \xb1*$ to produce an accurate prediction of the osmotic pressure in articular cartilage as it varies monotonically with increasing bath concentration. While the framework presented in this study could be applied to divalent (e.g., $Ca\u2009Cl2$) or trivalent (e.g., $Al\u2009Cl3$) cations, it should be emphasized that experiments would be needed to validate the osmotic response predicted from our model or any other. In addition, although the general multiphasic framework [12,73] may be used for neutral solutes and can account for their osmotic pressure with a suitable constitutive relation for the osmotic coefficient, the model formulated in Sec. 2.1.2 considers Donnan osmotic pressure generated by charged solutes alone and thus cannot be used for neutral solutes.

In summary, this study has presented direct measurements of the osmotic swelling pressure in bovine and human articular cartilage in buffered saline at various concentrations and described the osmotic pressure using a nonideal Donnan model with a constitutive relation for the osmotic coefficient which depends on a single parameter *ξ*. Finite element curve-fitting of the experimental transient stress-relaxation responses of these cartilage samples, using the calculated osmotic coefficient, exhibited a concentration-dependent aggregate modulus *H _{A}*, consistent with prior literature findings [15,23]. In contrast to the Poisson–Boltzmann theoretical interpretation proposed in earlier studies [15,18,23], which attributed the osmotic pressure as well as the concentration-dependent modulus to electrostatic interactions mediated by charged ions in the bathing solution [15], our modeling approach explicitly attributes the Donnan osmotic pressure to the osmolarity difference between the tissue and its surrounding bath resulting from ion transport in the presence of a charged solid matrix, as embodied in the triphasic theory [10,12]. Then, additional charge repulsion effects are accounted for by allowing the solid matrix material properties to vary with strain and ion concentrations, an approach consistent with the principle of equipresence in mechanics, which states that all functions of state must depend on all state variables unless explicitly precluded by the axiom of entropy inequality. We agree with the hypothesis that the Poisson–Boltzmann model is appropriate for describing the physics of the concentration dependence of

*H*in relation to charge-to-charge repulsion forces; however, the triphasic theory requires a constitutive relation of the form given in Eq. (3.1) to account for this effect. The results of this study expand the available experimental data on measurements of the osmotic pressure in bovine and human cartilage and provide a modeling approach that can capture the full complexity of this tissue in response to mechanical and chemical (osmotic) loading. This framework provides a valuable tool for computationally modeling complex electrokinetic phenomena in cartilage, thereby facilitating future work in modeling cartilage damage and fatigue processes which are often associated with tissue swelling. More broadly, models of the mechanical interactions between ECM constituents in other biological tissues may benefit from the nonideal osmotic behavior and models presented here.

_{A}## Acknowledgment

The authors would also like to acknowledge help from Glendon Zimmerman in developing software to process and analyze the thousands of curve-fits produced by finite element analysis, and Jalisa Zimmerman for fruitful discussions on statistical analysis and visualization.

## Funding Data

National Institute of General Medical Sciences (Grant No. GM083925; Funder ID: 10.13039/100000002).

National Science Foundation Graduate Research Fellowship Program (DGE-11-44155; Funder ID: 10.13039/100000001).

## Footnotes

In a linear theory, note that $\sigma e=HA\xb7\epsilon $ where *H _{A}* is the aggregate modulus of the tissue, which may depend on the ion concentrations.

In cartilage, the fixed charge density of proteoglycans is negative. Here, we use a positive value for *c ^{F}* and adjust the sign of our formulas to account for the negative charge.

Assuming the external bath electric potential is zero, the partition coefficient in the bath is unity, and the activity coefficients in the bath are both unity.