## Abstract

This study presents a damage mechanics framework that employs observable state variables to describe damage in isotropic or anisotropic fibrous tissues. In this mixture theory framework, damage is tracked by the mass fraction of bonds that have broken. Anisotropic damage is subsumed in the assumption that multiple bond species may coexist in a material, each having its own damage behaviour. This approach recovers the classical damage mechanics formulation for isotropic materials, but does not appeal to a tensorial damage measure for anisotropic materials. In contrast with the classical approach, the use of observable state variables for damage allows direct comparison of model predictions to experimental damage measures, such as biochemical assays or Raman spectroscopy. Investigations of damage in discrete fibre distributions demonstrate that the resilience to damage increases with the number of fibre bundles; idealizing fibrous tissues using continuous fibre distribution models precludes the modelling of damage. This damage framework was used to test and validate the hypothesis that growth of cartilage constructs can lead to damage of the synthesized collagen matrix due to excessive swelling caused by synthesized glycosaminoglycans. Therefore, alternative strategies must be implemented in tissue engineering studies to prevent collagen damage during the growth process.

## 1. Introduction

All biological tissues sustain some amount of fatigue or injury damage over time, even under normal physiological conditions. Most tissues also undergo continuous repair, though the corresponding rate of remodelling may not always keep up with the rate of damage, depending on various factors, including ageing and extent of injury. The ability to model the progression of mechanical damage in biological tissues may be valuable for better understanding the balance between damage and repair, and for elucidating the influence of pathomechanics of mechanically mediated progressive degenerative diseases on tissue composition and structure. For example, in articular cartilage, mechanical damage induced by loading may be the initiating factor in progressive osteoarthritic degeneration.

In tissue mechanics, theoretical models of damage have been applied most extensively for the analysis of bone [1–3] but also for soft tissues such as tendon [4] and cartilage [5]. These models are generally based on the continuum damage mechanics framework originally developed by Kachanov [6] and Rabotnov [7], and subsequently refined by LeMaitre [8,9] and Chaboche [10], Simo & Ju [11], and others. As summarized in [11], this framework for continuum damage mechanics is based on the thermodynamics of irreversible processes and internal state variable theory [12], where damage is modelled as an internal or hidden (i.e. non-observable) state variable.

Therefore, modern measurement techniques that can directly detect chemical bond modifications and disruptions in biological tissues, such as biochemical measurements for assessing cross-link density in collagenous tissues [13], immunohistochemistry for assessing mechanically or enzymatically degraded collagens [14–19], and Raman spectroscopy for detecting alterations in chemical bonds [20,21], may not inform the traditional damage framework directly.

The objective of this study was to demonstrate that this classical continuum damage mechanics framework may be reformulated using only observable scalar state variables for modelling damage in isotropic and anisotropic materials. These state variables represent the density of intact or broken bonds in the solid constituents of a tissue. Different bond species may exist in a tissue, each of which may regulate different components of the strain energy density function governing the anisotropic material response. In this framework, the evolution of bond densities with time is governed by the axiom of mass balance, with a mass supply term representing the exchange of mass between intact and broken bonds. The reaction representing the breaking of bonds is regulated by a failure criterion measure, such as strain energy density or effective stress, whose magnitude determines the probability of bonds breaking, based on a given probability density function (p.d.f.). The evolution of tissue ultrastructure with increasing bond damage may be taken into account by distinguishing the bond densities of fibril bundles oriented along different directions and allowing them to evolve independently based on each fibril bundle's state of strain. This approach differs from the alternative method of allowing the texture tensors of anisotropic materials from evolving with damage.

We show that this reactive damage framework reproduces the classical damage mechanics formulation for isotropic materials, with the classical damage variable *D* actually representing the fraction of bonds that are broken. Using a bond-breaking framework for damage mechanics demonstrates that damage may still be represented by a scalar variable *D* even in the case of anisotropic materials. Because of its relevance to soft tissue mechanics, the damage response of fibrous materials is examined. First, the damage response of a single fibre bundle is analysed for various choices of p.d.f.'s, under homogeneous and inhomogeneous strain states. Then, the damage mechanics of continuous fibrous distributions are examined, where bundles oriented along different directions are associated with direction-dependent damage variables.

Finally, this damage framework is used to examine an intriguing observation from our cartilage tissue engineering studies, regarding the temporal evolution of the composition and material properties of tissue constructs. Experiments reported here have shown that, in long-term culture, collagen and proteoglycan content grow continuously over time, even though the construct's compressive Young's modulus reaches a plateau. We use our damage framework to examine the hypothesis that collagen fibrils synthesized early in culture become progressively more damaged, as the increasing proteoglycan content causes greater Donnan osmotic swelling, thus shifting the structural integrity of the collagenous matrix to younger, less damaged, generations of synthesized collagens.

## 2. Theoretical formulation

Biological tissues often consist of heterogeneous mixtures of constituents. For example, in soft connective tissues, the solid matrix may consist of a variety of collagen types, elastin and other macromolecular species such as proteoglycans and hyaluronan. Each specific substance may exhibit a variety of bond species, which may get damaged under different conditions. Therefore, the theoretical formulation first addresses the damage response of a single-bond species, followed by a generalization to multiple bond species. Furthermore, this presentation only considers putative strong bonds that break permanently, thus producing damage, in contrast with putative weak bonds that break and reform spontaneously, thus producing viscoelasticity [22,23]. Lastly, damaged bonds in biological tissues may ultimately be reincorporated (rebound) due to a repair response. Such phenomena are not featured here, although such behaviours are particularly well accommodated within the theoretical framework. The framework for this theoretical formulation is the theory of mixtures [24], as applied specifically to constrained mixtures of solid constituents [25,26].

### 2.1. Single intact bond species

#### 2.1.1. Bond energy

Consider an intact bond species *b* whose specific bond energy (bond energy per mass) is given by , where **F** is the deformation gradient of all solid constituents. Let be the referential mass density (mass per volume in the reference, or undeformed, configuration) of all *intact* bonds of that species. Then, represents the referential strain energy density (strain energy density per volume in the reference configuration) when all bonds are intact. When some of the bonds have failed, the current referential density of remaining intact bonds is , such that , and the current referential strain energy density is given by
2.1
Thus, the strain energy density of the damaged bond species has been scaled from by the fractional ratio of remaining intact bonds to all the original intact bonds, . In this notation, the subscript 0 implies that all bonds are intact for that variable, whereas the subscript ‘r’ implies that only some of the bonds remain intact; for both subscripts, the corresponding variable is normalized by the volume of the material in the reference configuration.

Now consider that bonds in that species have a probability of failure determined by a failure criterion measure , where the definition of follows the classical framework for damage [27]. In this study, it is proposed that the probability that intact bonds of species *b* will fail at a value of is given by a p.d.f., . The choice of p.d.f. is based on a constitutive assumption, with the requirement that the p.d.f. be defined on the appropriate interval for , e.g. the semi-infinite interval when is a positive semi-definite function of **F**. The corresponding cumulative distribution function (c.d.f.) of is .

Let the failure measure be given by at a particular state of deformation **F**(*s*), and let the maximum value of over the history up until the current time *t* be given by [27]
2.2
Then, based on the formulation proposed in this study,
2.3
represents the fraction of all bonds that have failed up until the current time *t*, leaving
2.4
as the remaining fraction of intact bonds. Thus, equation (2.1) may be rewritten as
2.5

#### 2.1.2. Thermodynamics

To examine the thermodynamics of this process, consider that damage represents a reaction that transforms intact bonds into broken bonds consisting of two or more fragment species *c*,
2.6

The evolving referential density of broken fragments is , with prior to any damage, and when all bonds *b* have broken. In this reactive damage framework, it is assumed that the specific strain energy of broken bond fragments is zero, ; therefore, the net free energy of the mixture of *b* and *c* species is simply
2.7
which is given in equation (2.5). (This approach precludes the modelling of plasticity, where bonds may break and reform in a loaded state.) The chemical potential of each species *a* = *b*,*c* may be evaluated from the standard relation [26],
2.8
which shows that the chemical potential is equal to the bond specific free energy, with and . This equivalence between chemical potential and specific free energy occurs because of the simplifying constitutive assumption that *ψ ^{b}* is not a function of any of the bond densities . (In physical chemistry of solutions, this same assumption produces

*ideal*solutions; we may thus describe our bond mixture model as an ideal mixture.) The rate at which the intact bond density changes over time is governed by the axiom of mass balance for that constituent [26], 2.9 where is the material time derivative of and is the referential mass density supply to

*b*from reactions with all other constituents in the mixture, such as the bond-breaking reaction of equation (2.6); a constitutive relation must be provided for . Similar expressions may be given for the broken fragments

*c*and, according to the axiom of mass balance for the mixture, . As shown previously in the context of mixture theory [26], a reaction can proceed spontaneously according to the Clausius–Duhem inequality if and only if 2.10

As for all *c* species, this inequality constraint reduces to , which further implies that (thus, ) as *ψ ^{b}* is a positive semi-definite function of the strain. Thus, the thermodynamic constraint simply states that intact bond species may break spontaneously with loading (as their mass supply must be negative) but may not reform spontaneously in this damage framework.

#### 2.1.3. Damage criterion

As described in the prior damage mechanics literature [27], a damage criterion needs to be defined such as
2.11
where *φ ^{b}* represents a damage surface whose tensorial normal is
2.12

The damage criterion of equation (2.11) presents the following alternatives: *φ ^{b}* < 0, if the current state of deformation is below the current damage threshold, or

*φ*= 0 and the deformation is receding from, or tangent to, the damage (), or increasing the damage (). Substituting equation (2.4) into the mass balance of equation (2.9) produces the damage evolution constitutive relation, 2.13 showing that the constitutive relation for is thus predicated on the shape of the c.d.f. . Based on this relation, the thermodynamic requirement, , and the damage progression conditions,

^{b}*φ*= 0 and (implying at ), combine to produce the thermodynamic constraint 2.14 on the c.d.f. Therefore, the entropy inequality is satisfied when is a monotonically increasing function over the domain of .

^{b}#### 2.1.4. Stress

For a single bond species *b*, the state of stress may be evaluated from using the usual relations of hyperelasticity, recalling that and **F** are independent state variables [26]. Thus, for the Cauchy stress,
2.15
where *D ^{b}* is evaluated from equation (2.3) and
2.16

### 2.2. Failure criterion measure

The failure criterion measure of a particular bond species *b* may be based on strain invariants, strain energy density, stress invariants or any other suitable, experimentally validated, measure. Measures that are dependent on the material behaviour, such as strain energy density and stress, are based on the behaviour of intact bonds, thus for the strain energy density, and for the Cauchy stress. For example, Hayhurst [28] proposed to use a linear combination of the maximum principal stress , the hydrostatic stress and the octahedral shear stress (or von Mises effective stress) of , to model creep damage, thus
2.17
where *α* and *β* are material constants and
2.18

Simo [27] proposed an equivalent strain measure, 2.19 deemed appropriate for isotropic damage. It is also possible to select . For applications to biological tissue damage these criteria may be associated with various mechanisms; selection of a specific criterion is dependent on the tissue's constitutive behaviour and the applicability of different criteria must be guided by experimental data.

#### 2.2.1. Cumulative distribution functions

The choice of p.d.f.'s and corresponding c.d.f.'s is predicated on experimental observation of damage. For example, the lognormal distribution is a suitable p.d.f. when , with
2.20
where *μ*_{c} (same units as ) and *σ*_{c} (unitless) represent material constants, such that half of the bonds break at (figure 1).

The Weibull distribution is a similarly suitable function, with
2.21
where *μ*_{c} (same units as ) and *α*_{c} (unitless) are material constants, such that the fraction 1 – *e*^{−1} of bonds break at (figure 2).

With decreasing *σ*_{c} and increasing *α*_{c}, these two distributions exhibit opposite biases about the value of *μ*_{c}: in the lognormal distribution, damage progresses more rapidly with increasing when , and more slowly when ; the opposite trend is found in the Weibull distribution. Alternatively, when modelling fracture (onset of complete damage at a single threshold value of ), the p.d.f. may be given by the Dirac delta function , with its c.d.f. given by the Heaviside unit step function ,
2.22
where *μ*_{c} is the critical value of at which fracture occurs.

### 2.3. Multiple intact bond species

This reactive damage framework may be generalized to any number of intact bond species *b*, each associated with its own strain energy density , such that the net free energy density of the mixture is
2.23
and the mixture Cauchy stress is
2.24
For each of these bond species, the damage may be given by a different c.d.f. and a different functional form of the failure criterion measure.

This approach can accommodate anisotropic materials and anisotropic damage, by allowing damage to evolve differently for various terms of the strain energy density function *Ψ*_{r}. For example, in a transversely isotropic material whose strain energy density depends on five strain invariants, damage may evolve according to
2.25
This form of splitting damage among various terms is only an illustration (a constitutive assumption) and does not represent a general form (which does not exist). As the molecular structure of a material permits any number of bond species, constitutive relations for damage may include any number of bond-breaking reactions in response to loading.

### 2.4. Comparison with classical damage mechanics

In isotropic materials, damage is represented by a scalar state variable *D*, where *D* = 0 in the undamaged state and *D* = 1 in the completely damaged state. More generally, for anisotropic materials, a tensorial state variable is needed. For example, in infinitesimal strain theory, a fourth-order tensorial function of the fourth-order tensorial damage is used to attenuate the intact (or virgin) stress according to
2.26
where represents the stress at the current state of strain, in the absence of damage, and **T** is the actual stress [29,30]. For isotropic materials, , where is the fourth-order identity tensor and **I** is the second-order identity tensor, and , so that , consistent with our formulation in equation (2.15), though derived using a different framework. A constitutive relation is also needed to model the evolution of damage, such that the damage time rate of change is given by
2.27

This theoretical framework requires the formulation of fourth-order tensor relations for and in the general case of anisotropic materials [29,30]. Furthermore, as an internal variable, is strictly non-observable [12]; therefore, it cannot be measured experimentally.

In our mixture framework where damage is tracked by the mass fraction of bonds that have broken, anisotropic damage is subsumed in the assumption that multiple bond species may coexist in a material, each having its own damage behaviour, as covered in the previous section and summarized in equations (2.23) and (2.24). The mixture approach employs bond kinetics to formulate the constitutive relation for in the mass balance relation of equation (2.9), whereas the damage evolution relation in equation (2.27) is not based on any of the classical balance axioms because it employs non-observable variables.

## 3. Fibrous solid matrix

### 3.1. Unidirectional fibre bundle

Consider a fibrous biological soft tissue whose ground substance is modelled as a compressible neo-Hookean solid whose free energy density is [31]
3.1
where and are material constants equivalent to the Lamé constants in the limit of infinitesimal strains, and , where **C** is the right Cauchy–Green tensor. Consider that the fibrous matrix consists of a unidirectional fibre bundle with a single bond species, whose free energy density is
3.2
where is a material constant related to the fibre stiffness, *β*_{f} is a unitless exponent (*β*_{f} ≥ 2), and is the square of the stretch ratio along the referential fibre direction **n**. The Heaviside unit step function *H*(·) ensures that the fibres contribute to the tissue response only when they are in tension (*I _{n}* > 1).

For illustrative purposes, consider that only the fibres' bonds may get damaged in this tissue, such that the tissue free energy density is
3.3
where *D* is the damage variable for the single-bond species of the fibres. Let the failure criterion measure be the fibre free energy density, , and let the damage c.d.f. either be the lognormal or the Weibull distribution. Finite-element analyses are performed using a custom implementation of this damage model in the open source finite-element program FEBio [32] (www.febio.org).

In the first analysis, consider that the fibrous tissue is subjected to a homogeneous state of uniaxial tensile stress with traction-free lateral boundaries, where the axial stretch ratio is increased from 1 to 2. The material properties for the ground substance are *μ*_{g} = 0.05 megapascals (MPa), *λ*_{g} = 0.075 MPa; and for the matrix, *ξ*_{f} = 1 MPa and *β*_{f} = 2. These properties are representative of a fibrous tissue such as articular cartilage [33–36]. The fibres are oriented along the loading direction. In the absence of damage, the value of *Ψ*_{r} peaks at 4.5 mJ mm^{−3} when the stretch ratio is 2, with a small contribution from the ground substance, ; thus, almost all the strain energy density is stored in the fibres.

When damage is modelled using a lognormal c.d.f., with *μ*_{c} *=* 1 mJ mm^{−3} and *σ*_{c} = 1/2, 1/4 or 1/8, the free energy density *Ψ*_{r} begins to deviate from the intact response at values below *μ*_{c} (figure 3). Smaller values of *σ*_{c} exhibit a later onset and quicker propagation of damage than larger values of this parameter, as evidenced in the profiles of *Ψ*_{r} and the axial stress. As the ground substance was not assumed to get damaged in this example, a small residual strain energy density and axial stress are found to persist even after complete failure of all fibres.

A similar response is observed when damage is modelled with a Weibull c.d.f., with *μ*_{c} = 1 mJ mm^{−3} and *α*_{c} = 2.5, 5 or 10, where larger values of *α*_{c} exhibit a later onset and quicker propagation of damage than smaller values of this parameter (figure 4).

### 3.2. Multiple fibre bundles

Many fibrous biological soft tissues exhibit a continuous fibre distribution density about one or more predominant fibre directions, including articular cartilage, arterial wall and skin [37–39]. As the tissue is loaded, fibres get recruited along the loading direction, producing a load response that is significantly dependent on the initial fibre distribution and recruitment [40]. To examine the damage response of such fibre distributions, first consider the example of a tissue with three coplanar fibre bundles, oriented at *θ _{b}* relative to the loading direction (

*b*= 1, 2, 3). The free energy density of this tissue is given by 3.4 where is given by equation (3.1) and has the form given in equation (3.2), with distinct values of the fibre properties for each bundle. Similarly, consider that the damage model for each fibre bundle follows the Weibull c.d.f. with , such that the damage variable

*D*is similarly associated with distinct values of the Weibull parameters .

^{b}As in the previous analysis, the fibrous tissue is subjected to a homogeneous state of uniaxial tensile stress with traction-free lateral boundaries, where the axial stretch ratio is increased from 1 to 2. Three cases are examined: (a) the three fibre bundles are all oriented along the same direction *θ _{b}* = 0 relative to the loading direction; they share the same fibre properties and ; however, the Weibull parameters differ such that and , 1 and 1.5 mJ mm

^{−3}for

*b*= 1, 2 and 3, respectively. (b) The fibre bundles are oriented at

*θ*= 0 and ±

_{b}*π*/6 whereas all three fibre bundles share the same properties, , , , . (c) The first bundle is oriented along

*θ*

_{1}= 0, with and and the other two (

*b*= 2, 3) are oriented along

*θ*

_{b}*=*±

*π*/6, with and ; other properties are the same for all three bundles, , .

The axial stress–stretch responses for these three cases exhibit distinct behaviours (figure 5*a*). In all cases, as one fibre bundle gets damaged, fibres from the remaining bundles sustain greater stress, showing progressive rises and falls in the stress–stretch profile. For case (a), three consecutive failure events may be noted in the response, indicating that the three aligned bundles contribute distinctly at different levels of the prescribed axial deformation; the peak stress reached with increasing damage in the second bundle is higher than that in the first bundle. For cases (b) and (c), only two consecutive failure events are observed, as two of the fibre bundles (*θ _{b}*

*=*±

*π*/6) are oriented symmetrically about the third bundle which coincides with the loading direction; thus, damage progresses identically in these symmetric bundles, explaining the presence of only two peaks. In (b), the second peak is higher than the first peak, whereas in (c) the opposite occurs. An average measure of the damage in all three fibres is also shown, evaluated as with

*N*= 3 (figure 5

*b*). The plateaus observed in the plots of versus axial stretch ratio indicate the occurrence of failure of fibre bundles; the subsequent rise in represents the increasing damage in the remaining fibre bundles being recruited, until they fail as well.

### 3.3. Discrete fibre distributions

The results from the above analysis of discrete fibre bundles suggest that a continuous fibre distribution may similarly require making constitutive assumptions about the elastic and damage properties along various directions. Normally, for continuous fibre distributions, variations in elastic properties with direction are taken into account using a fibre density distribution function *R*(**n**), where **n** represents the fibre direction in the reference configuration [39,40]. The function *R*(**n**) typically scales the free energy density and stress responses such that
3.5
where the integration is performed over the unit sphere, for three-dimensional distributions, or the unit circle, for two-dimensional distributions. The fibre density distribution function must satisfy . In the integrand, represents the strain energy density of the fibre bundle oriented along **n**. For example, may have the form given in equation (3.2).

To account for distinct damage responses along each fibre bundle direction **n**, let the free energy density of the tissue be given by
3.6
which represents a generalization of equation (3.4) to continuous fibre distributions. The function *D*(**n**) follows a c.d.f. whose damage properties may vary with **n**, based on experimentally validated constitutive assumptions. For example, for any of the c.d.f.'s presented in equations (2.20)–(2.22), it may be assumed that the parameter varies with the direction **n** according to , where is a constant. As *R*(**n**) scales the intact fibre response in a continuous fibre distribution, this assumption for effectively proposes that directions along which fibres exhibit lower (higher) stress responses at a given *I _{n}* also get damaged at lower (higher) values of the failure criterion measure .

In practice, the integral of equation (3.6) may be approximated by a summation,
3.7
where the unit sphere (for three-dimensional distributions) or unit circle (for two-dimensional distributions) is discretized into *N* finite elements, each with area Δ*A ^{i}*, whose centroid is located along

**n**

*, with . A variety of methods are available for discretizing the unit sphere into finite elements with nearly uniform Δ*

^{i}*A*values across the surface [41,42], including a geodesic dome representation [43]. We may refer to the discretized form in equation (3.7) as a discrete fibre distribution (d.f.d.). In the absence of damage (

^{i}*D*(

**n**

*) = 0), this summation converges to the continuous fibre response as*

^{i}*N*→ ∞.

A common choice for *R*(**n**) in biological soft tissues is the *π*-periodic von Mises distribution [44]. For three-dimensional fibre distributions with a principal orientation along the unit vector **m**,
3.8
where *β* > 0 is the fibre concentration parameter and is the imaginary error function. In the limit as *β* *→* 0 this distribution becomes isotropic, with . For two-dimensional fibre distributions with a principal orientation along **m**, the equivalent expression is
3.9
where and **n** lies in the plane of the fibre distribution; *I*_{0} is the modified Bessel function of the first kind of order 0. In the limit as *β* *→* 0, this distribution becomes transversely isotropic, with .

To illustrate the damage response of continuous three-dimensional fibre distributions, consider the model of equation (3.6) where is given in equation (3.1), has the form given in equation (3.2), and *R*(**n**) is given by equation (3.8). Let *D*(**n**) follow the Weibull c.d.f. of equation (2.21), with and , where is a constant. As in the previous analyses, the fibrous tissue is subjected to a homogeneous state of uniaxial tensile stress with traction-free lateral boundaries, where the axial stretch ratio is increased from 1 to 10. The material properties for the ground substance are *μ*_{g} = 0.05 MPa, *λ*_{g} = 0.075 MPa; and for the fibres, *ξ*_{f} = 1 MPa and *β*_{f} = 2. For the Weibull c.d.f., and *α*_{c} = 10. Two different fibre concentrations are considered in these examples: isotropic (*β* = 0) and transversely isotropic with a relatively high concentration (*β* = 5), with **m** coinciding with the axial loading direction.

Using the discretization of equation (3.7), the number *N* of fibre bundles needed to reproduce the conditions of a continuous fibre distribution is first investigated for the case *β* = 0, using *N* = 196, 596, 1196 and 1796. Results show that the axial stress versus stretch response is nearly identical for all *N* values up to a stretch ratio of approximately 2, when no significant damage has yet occurred, as also suggested by the response of an intact (i.e. non-damageable) continuous fibre distribution (figure 6*a*). This outcome indicates that the stress calculation in the absence of damage has converged with increasing *N*, producing the solution for a continuous fibre distribution. Beyond this stretch ratio, however, the onset of damage varies with *N*, such that the peak value of the axial stress (a measure of damage resilience) increases with increasing *N*. A plot of the average damage in the fibres shows that increasing *N* decreases the average damage per fibre at a given stretch ratio (figure 6*b*).

A similar behaviour is observed in the case of a concentrated transversely isotropic fibre distribution (*β* = 5), with the main distinctions being that the stress–stretch response exhibits a significantly higher effective modulus relative to the isotropic fibre distribution, and the progression of damage occurs more rapidly with increasing stretch, after the peak axial stress has been achieved (figure 7).

These results demonstrate that the damage response is sensitive to the number of fibre bundles; increasing the number of fibre bundles increases the resilience of the fibrous tissue to damage. A plot of the peak axial stress achieved for each *N* indicates that this relation is nearly linear, for both isotropic and transversely isotropic fibre distributions (figure 8). As the damage resilience increases as , a truly continuous fibre distribution (where *N* → ∞) would exhibit damage only as the stress goes to infinity. We conclude that seeking a converged numerical solution for the damage response as *N* → ∞ is nonsensical: damage is in fact critically dependent on the number of fibre bundles. In these two illustrations, fibres oriented along different directions increasingly get recruited toward the loading direction with increasing stretch and damage, providing further resilience.

## 4. Cartilage tissue engineering

In the experimental component of this study, we examine the topic of collagen damage in a growing engineered cartilage tissue construct. The concept that growth may cause damage is perhaps counterintuitive, but this hypothesis was motivated by experimental observations from our tissue engineering studies, as described here.

Cartilage tissue engineering combines cells and a biocompatible scaffold to synthesize tissue constructs for the purpose of repairing defects in osteoarthritic articular layers. When cultured in a chemically defined medium, constructs synthesize native extracellular matrix proteins, principally glycosaminoglycans (GAG) and collagen. Sulfated GAG imparts a negative fixed-charge density to the tissue, causing it to swell due to the resulting Donnan osmotic pressure; conversely, a cross-linked collagen network functions to resist this osmotic swelling [45]. Together, the collagen network and GAG give rise to the robust mechanical properties observed in native cartilage.

In our model system, we have used 2% agarose hydrogel as the scaffold and primary juvenile bovine chondrocytes as the cell source [46]. In recent years, we have adopted a chemically defined (serum-free) chondrogenic medium, supplemented with 10 ng ml^{−1} transforming growth factor (TGF)-*β*3 for the first 14 days in culture [47]; this temporary supplementation of growth factor has been shown to double the rates of GAG and collagen synthesis after day 14 [48,49]. This methodology has consistently produced GAG levels that meet or exceed those of native juvenile bovine cartilage, producing equilibrium compressive Young's moduli *E*_{Y} that similarly meet or exceed levels in native tissue [47]. For example, in our recent study using a chondrocyte seeding density of 30 × 10^{6} cells ml^{−1}, GAG content was 5.7 ± 1.6% of wet weight (%ww) and Young's modulus was 483 ± 65 kPa at day 45 [49]. In comparison, native GAG content in juvenile bovine cartilage is 2.8 ± 0.2%ww and *E*_{Y} = 400 ± 200 kPa [50]. By contrast, collagen levels in our tissue constructs have fallen consistently below native levels [47,49]; in the same tissue engineering study quoted above, collagen content was 2.4 ± 0.8%ww [49], whereas the native level is 12.3 ± 0.4%ww [50]. The relatively low level of collagen in our constructs was reflected in their low tensile equilibrium modulus (*E*_{T} = 703 ± 315 kPa) when compared with native tissue levels (*E*_{T} = 17 ± 5 MPa [51]).

To address this deficiency in collagen content, we previously investigated the effect of digesting GAG using chondroitinase-ABC (CABC) at intermediate time points in culture, to let collagen synthesis catch up with GAG synthesis [52,53]. These results produced higher collagen content than undigested controls, but did not reach native levels.

In this study, we opted to increase the cell seeding density to 120 × 10^{6} cells ml^{−1}, to match native levels in immature tissue (100–300 × 10^{6} cells ml^{−1} [54]). We also increased media volume per construct to account for the nutrient needs of such elevated cell numbers, consistent with our recent findings [55]. Finally, we extended the culture duration to 105 days to provide more time for matrix synthesis, and included a CABC treatment group with the digestion occurring from day 14 to 16 using 0.15 U ml^{−1} as per our previous protocol [53]. Cylindrical constructs had an initial diameter of 3 mm and height of 2.3 mm. GAG, collagen and cell content, and equilibrium compressive Young's modulus, were assessed using standard protocols described in our previous studies [55]. Two-way analysis of variance (ANOVA) (*α* = 0.05) was used to detect statistical differences (*p* ≤ 0.05) for the factors of treatment (control and CABC) and time (14/16, 35, 56, 77 and 105 days for control/CABC), with *n* = 4–6 samples per group; Tukey–Kramer correction was used for *post hoc* testing of the means.

One of the most striking outcomes of this experimental study was the considerable tissue construct volume growth or expansion from the size of the day-0 construct (figure 9). The control group specimens saw a nearly ninefold increase in their wet weight (and volume) relative to day 0, whereas CABC-treated specimens increased more than fivefold (figure 10). Because of this very large swelling ratio *J*, specifically defined as construct wet weight at a given time point normalized to the day-0 construct wet weight, GAG and collagen compositional measures are reported using normalization by current day wet weight (final day in culture) to produce conventional fractional measures (%ww, figure 11), and normalization by day 0 wet weight to assess total matrix synthesis (%D0ww, figure 12).

Using conventional measures, GAG content rose significantly above day 14/16 levels to nearly 7%ww by day 56 (*p* < 10^{−4}), then dropped slightly by day 105 (figure 11*a*); differences were observed between control and CABC treatments only at day 14/16 (*p* < 10^{−3}). Collagen content rose significantly above day 14/16 levels to nearly 2%ww by day 56, then remained unchanged until day 105 (figure 11*b*); no differences were observed between control and CABC treatments (*p* = 0.29).

Using day 0 normalization, GAG content increased nearly linearly with time (*p* < 10^{−4}), reaching 54%D0ww in the control group and 31%D0ww in the CABC group at day 105 (figure 12*a*). Collagen content similarly rose nearly linearly with time (*p* < 10^{−4}) to 15%D0ww in the control group and 10%D0ww in the CABC group at day 105 (figure 12*b*). From day 56 onward, GAG (*p* < 10^{−4}) and collagen (*p* < 10^{−4}) content in the CABC group were smaller than control.

Total cell content per construct did not exhibit a monotonic trend, though significant differences were found with time in culture (*p* < 10^{−4}); the CABC group exhibited consistently lower cell content than control (*p* < 10^{−4}) (figure 13). Fractional fluid content in constructs decreased over time until day 35 (*p* < 10^{−4}), then remained nearly constant until day 105 (figure 14); at most time points, CABC constructs had lower fluid content than controls (*p* < 10^{−4}). Finally, the compressive equilibrium Young's modulus rose significantly above day 14/16 levels until day 35, then remained statistically constant until day 105 (figure 15). At two time points, *E*_{Y} was higher in the CABC group compared with control (*p* < 0.01).

## 5. Damage in growing tissue constructs

On first pass, results of this experimental study were substantially consistent with our past results [49,52], producing comparable GAG and collagen %ww levels and compressive moduli, despite the fact that the cell seeding density was considerably higher. No differences were found in the conventional GAG and collagen measures between control and CABC treatments, though *E*_{Y} was significantly higher in the CABC group.

However, once we accounted for the swelling ratio, our interpretation of these results was considerably refined: using a cell seeding density comparable with native juvenile tissue, it was in fact possible to synthesize native levels of collagen in our tissue constructs (15%D0ww in the control group at day 105), indicating that the primary chondrocytes used in this culture system had not lost their ability to produce sufficient collagen. However, the supra-physiological levels of GAG synthesized over the same time period (54%D0ww, compared with the native GAG content of 2.8%ww) produced so much swelling, resulting from the excessive GAG deposition-induced swelling pressure, that the collagen content was diluted (2%ww). It was experimentally evident that the swelling was caused primarily by the GAG (consistent with Donnan swelling theory [56]), as lower GAG %D0ww levels in the CABC group (figure 12*a*) produced a lower swelling ratio (figure 10).

The next intriguing outcome was that *E*_{Y} was higher in the CABC treatment than control, even though %ww levels of GAG and collagen were not different, and %D0ww levels were in fact smaller in the CABC group. As GAG %ww usually correlates strongly with *E*_{Y} in tissue constructs [57], this surprising outcome led us to the hypothesis that the large swelling ratios observed in this study could in fact damage the collagen network, whose integrity is critical to prevent uncontrolled swelling (as the agarose hydrogel breaks apart early in the growth process based on visual inspection).

We hypothesized the following specific mechanism: GAG and collagen are synthesized continuously during the growth process, each at nearly constant rates. The charged GAG produces a swelling pressure that expands the tissue construct; this expansion is primarily resisted by the collagen network that was synthesized at the start of culture until the current time. Collagen synthesized at earlier times points is then subjected to more swelling than collagen synthesized most recently, under the assumption that collagen is incorporated into the network in a stress-free (unswollen) state. Therefore, as time progresses, earlier generations of collagen are subjected to considerable stretching (up to ninefold swelling for the earliest collagen generation in this study, figure 10), causing them to become damaged. Only the most recent collagen generations remain intact or nearly intact, providing some resistance to swelling while maintaining a regular construct shape (figure 9). At some time point during the culture (e.g. day 35 based on the data), the rate of matrix synthesis and damage equalizes, such that GAG and collagen %ww, and *E*_{Y}, all reach a plateau.

Based on our earlier study [49], nearly 90% of synthesized GAG and 85% of synthesized collagen is retained in the constructs throughout the culture period; therefore, in our hypothesized mechanism, we expect that damaged collagen mostly remains in the construct, possibly contributing to the structural integrity when newly synthesized generations cross-link with undamaged segments of those collagen fibrils.

Under this hypothesized mechanism, the CABC-treated constructs are subjected to less damage in their earlier generations because of the reduced swelling resulting from the GAG digestion. Therefore, the collagen network in the CABC group has a higher tensile stiffness, which effectively translates to greater resistance to compression from a swollen state as shown in our earlier studies on native cartilage [36].

To validate this hypothesized mechanistic interpretation of our experimental observations, we performed a theoretical analysis that captured its salient points. Cartilage constructs were modelled as a mixture of collagen, charged GAG, and an interstitial fluid containing water and monovalent salt counter-ions (Na^{+} and Cl^{−}). In this case, the mixture stress is
5.1
where *p* is the interstitial fluid pressure, **I** is the identity tensor, and **T*** ^{e}* is the collagen matrix stress. Under steady-state conditions, when interstitial fluid flow and ion diffusion have subsided, the interstitial fluid pressure is given by the Donnan osmotic pressure relative to ambient bath conditions,
5.2
where

*R*is the universal gas constant,

*θ*is the absolute temperature,

*Φ*is the osmotic coefficient describing the deviation from ideal Donnan law (

*Φ*= 1 for ideal conditions),

*c**is the external NaCl bath concentration and

*c*

^{F}is the fixed-charge density arising from the GAG. The parameter

*c*

^{F}was calculated directly from the measured GAG %D0ww using 5.3 where

*z*

^{GAG}is the charge number and

*M*

^{GAG}is the molar mass of GAG (

*z*

^{GAG}= −2 and M

^{GAG}= 513 g mol

^{−1}, recognizing that GAG content was measured using a dimethylmethylene blue assay with chondroitin sulfate as the standard);

*ρ*

_{T}is the construct mass density (taken to be

*ρ*

_{T}= 1 g ml

^{−1}) and

*φ*is the measured fractional fluid content. The only adjustable parameter in this relation is thus

^{w}*Φ*.

The collagen matrix stress **T*** ^{e}* was evaluated using the multigeneration growth model formulated in our earlier study [26]: each generation of collagen was assumed to be deposited in a stress-free reference configuration

**X**

*, coinciding with the current configuration at the start of the generation, , where*

^{γ}*χ*is the motion,

**X**is the reference configuration of the first generation (start of day 1 in culture), and

*t*is the time at the start of generation

^{γ}*γ*. The deformation gradient of each generation is thus and the collagen in generation

*γ*is subjected to strains calculated from

**F**

*. The construct swelling ratio is given by*

^{γ}*J*= det

**F**, where is the deformation gradient of the first generation. Though collagen synthesis is continuous, a computational scheme employs finite time increments and we chose to have 105 generations, corresponding to the start of each day in culture.

Each collagen generation *γ* was modelled using a fibre distribution with *N* = 210 fibre bundles corresponding to an isotropic distribution (*R*(**n*** ^{i}*) = 1/4

*π*in equation (3.7)), reflecting the random arrangement of collagen fibres evident in free-swelling tissue constructs. Each bundle was modelled using the constitutive relation of equation (3.2), with

*β*

_{f}= 2 and

*ξ*

_{f}left as a single adjustable parameter (same

*ξ*

_{f}value for all bundles and all generations). The lognormal c.d.f. of equation (2.20) was used to describe the probability of damage in each bundle of each generation, with adjustable parameters

*μ*

_{c}and

*σ*

_{c}(same values for all bundles). The failure criterion measure was the Simo criterion of the fibre bundle, equation (2.19), calculated for each collagen generation using

**F**

*for that generation.*

^{γ}This model was implemented in the FEBio finite-element environment to predict construct swelling and damage, and construct equilibrium compressive modulus under conditions equivalent to the experimental testing protocol (unconfined compression at 10% strain). Using FEBio predictions, the four adjustable material parameters, *Φ*, *ξ*_{f}, *μ*_{c} and *σ*_{c}, were fitted simultaneously to the mean experimental results for swelling ratio *J* versus time in culture (figure 10) and compressive modulus *E*_{Y} versus time in culture (figure 15), using control group data only (*Φ* = 0.5, *ξ*_{f} = 31 kPa, *μ*_{c} = 0.08(mJ mm^{−3})^{0.5} and *σ*_{c} = 1.7). These fits produced coefficients of determination and , respectively. To validate the model, these optimal values were then used to predict *J* and *E*_{Y} for the CABC group, producing coefficients of determination and , respectively (figures 10 and 15).

The successful curve-fits of control group data, and faithful theoretical predictions of CABC group data, provide a mechanistic validation of our hypothesis that growth is causing collagen damage in our engineered tissue constructs. The fact that the model could predict an initial rise in *E*_{Y} followed by a plateau (figure 15), whereas the swelling ratio *J* rose continuously (figure 10), indicates that the onset of damage with growth can explain the observed responses.

The agreement between theory and experimental datasets (control and CABC groups) is also evident in a scatter plot of *J* versus GAG %D0ww data, showing a strong correlation with the theoretical model (figure 16*a*). A scatter plot of *E*_{Y} versus GAG %ww from the control and CABC groups here and our prior study of [49] produces a moderate correlation (figure 16*b*).

## 6. Observable damage

Based on the general agreement between experimental results and our hypothesized damage-mechanism model of tissue growth, we investigated methods of measuring the damage directly and quantitatively. In 1997, Bank *et al.* [58] introduced a simplified measurement of collagen fibre denaturation using *α*-chymotrypsin, which reports the fraction of all collagen in that tissue that has been degraded (i.e. whose triple helix is unwound). They demonstrated that this assay could predict higher collagen degradation in cartilage from human osteoarthritic joints (10.4 ± 5.6%) versus normal tissue (2.1 ± 1.0%) suggesting that the technique would serve as an indicator of for collagen fibre damage.

We used this assay to assess damage in the control and CABC groups of the current study, as well as damage in native juvenile bovine cartilage from which we harvested chondrocytes for our tissue engineering studies. Native juvenile cartilage (*n* = 5 samples) exhibited 10.1 ± 0.9% denatured collagen, consistent with the approximately 14% value reported previously for young adult bovine cartilage [17]. Engineered cartilage constructs exhibited damage that increased over time, ranging from 38.9 ± 4.8% at day 14 to 53.9 ± 5.7% at day 105 for the control group, and from 22.4 ± 4.6% at day 16 to 63.8 ± 4.4% at day 105 for the CABC group (figure 17). A two-way ANOVA for the factors of time in culture and treatment showed significant differences across all factors (*p* < 10^{−4} for time, indicating a significant increase in damage over time; *p* = 0.002 for treatment, indicating that damage was lower in the CABC group versus control on average; and *p* < 10^{−4} for day × treatment, indicating that damage was lower in the CABC group only on day 14/16, figure 17). These experimental results confirmed that damage did occur in tissue constructs, well above the baseline level of denatured collagen found in native juvenile cartilage.

For comparison purposes, the fractional damage predicted from the model is also shown in figure 17). This damage was calculated by taking into account the multigenerational growth process: the average of over all *N* fibre bundles of generation *γ* was denoted by , which represents the damage in collagen generation *γ*; as each day produced a new generation, the fractional damage at the end of day *d* was evaluated from the cumulative damage until that day,
6.1

The trend of this theoretical versus time in culture was consistent with the biochemical assay results for both control and CABC groups, though the magnitudes did not match exactly (figure 17). This qualitative agreement was similarly supportive of our hypothesis that growth and swelling cause damage in the collagen matrix of engineered cartilage constructs. The fact that *D* did not agree exactly with the measured collagen degradation was not surprising, as the model did not propose an explicit constitutive relation between the fraction of collagen that is denatured and the fraction of molecular bonds that are broken.

## 7. Discussion

The objectives of this study were to formulate a damage mechanics framework that employs observable state variables, instead of internal hidden variables as in the prior literature, while recovering the classical formulation for damage in isotropic materials; to show that this new framework was applicable to anisotropic materials without appealing to a tensorial damage measure; to demonstrate its application to fibrous tissues with unidirectional or multi-directional fibre orientations and to use this framework to test the hypothesis that growth of cartilage constructs can lead to damage of the synthesized collagen matrix.

Most biological tissues exhibit viscoelasticity, which significantly hinders the ability to assess their damage exclusively from mechanical responses (except when the tissue exhibits gross failure). Similarly, damage assays (e.g. biochemistry or Raman spectroscopy) are not always easy to interpret on their own, as they may depend on complex interactions that are not fully understood. Using observable variables to track damage implies that predictions from damage models may be validated experimentally, and experimental damage measures may be interpreted mechanistically by validated models. This combination of theory and experiments can enhance our confidence when testing hypotheses, as demonstrated in the tissue engineering study presented here.

Although we have conducted cartilage tissue engineering studies for nearly two decades, we had not initially suspected that newly synthesized extracellular matrix could rapidly become damaged during the growth process. In retrospect, the first indication that excessive GAG swelling could cause damage to engineered cartilage was reported in our earlier study, where constructs were found to crack at their center [59], a one-time finding that we attributed to special circumstances. In this study, no cracks were observed in the constructs, but the temporal evolution of construct mechanical properties and swelling led us to hypothesize a more pervasive collagen damage process that we proceeded to investigate more quantitatively.

Had we only performed the biochemical assay (using *α*-chymotrypsin digestion [58]), our results would have remained somewhat ambiguous because this assay had not been previously used on engineered tissue constructs. Thus, the elevated levels of denatured collagen found here (figure 17) could have been attributed to the fact that newly synthesized cartilage has relatively high levels of immature collagen that may be responsible for these findings. (Indeed, it was somewhat surprising to us that healthy native juvenile and young adult bovine cartilage exhibit between 10% and 14% damage [17], whereas healthy adult human cartilage only exhibits 2% damage [58]). It was the combination of modelling and experiments that increased our confidence in the testing of our hypothesis.

The results of this study have significant implications for our ongoing cartilage tissue engineering efforts: increasing cell seeding density to levels comparable with neonatal cartilage can produce collagen synthesis that matches native levels; however, we must develop tissue engineering strategies that prevent damage of newly synthesized collagens. This may be achieved by either preventing excessive swelling during the growth process (e.g. suitably constraining the constructs without compromising access to nutrients), or reducing GAG synthesis in the early stages of tissue growth using alternative strategies to CABC (e.g. small interfering RNA to temporarily suppress aggrecan synthesis). These alternative strategies are currently under consideration.

## Data accessibility

Experimental data for engineered cartilage constructs reported in this manuscript are available for download from https://academiccommons.columbia.edu/catalog/ac%3A187860 and the model is available on the FEBio forum at http://mrlforums.sci.utah.edu/forums/.

## Competing interests

We declare we have no competing interests.

## Funding

Research reported in this publication was supported by the National Institute of Arthritis, Musculoskeletal and Skin Diseases and the National Institute of General Medical Sciences of the National Institutes of Health under award nos R01AR060361, R01AR043628 and R01GM083925.

## Disclaimer

The content is solely the responsibility of the author and does not necessarily represent the official views of the National Institutes of Health.

## Footnotes

One contribution of 19 to a theme issue ‘Integrated multiscale biomaterials experiment and modelling: towards function and pathology’.

- © 2015 The Author(s)