## Abstract

The structural integrity and the biomechanical characteristics of ligaments and tendons result from the interactions between collagenous and non-collagenous proteins (e.g. proteoglycans, PGs) in the extracellular matrix. In this paper, a dissipative theory of temporary interfibrillar bridges in the anisotropic network of collagen type I, embedded in a ground substance, is derived. The glycosaminoglycan chains of decorin are assumed to mediate interactions between fibrils, behaving as viscous structures that transmit deformations outside the collagen molecules. This approach takes into account the dissipative effects of the unfolding preceding fibrillar elongation, together with the slippage of entire fibrils and the strain-rate-dependent damage evolution of the interfibrillar bridges. Thermodynamic consistency is used to derive the constitutive equations, and the transition state theory is applied to model the rearranging properties of the interfibrillar bridges. The constitutive theory is applied to reproduce the hysteretic spectrum of the tissues, demonstrating how PGs determine damage evolution, softening and non-recoverable strains in their cyclic mechanical response. The theoretical predictions are compared with the experimental response of ligaments and tendons from referenced studies. The relevance of the proposed model in mechanobiology research is discussed, together with several applications from medical practice to bioengineering science.

## 1. Introduction

Ligaments and tendons are densely packed, regularly arranged, soft connective tissues having the biomechanical function both to transfer forces and to guide joint movements. Ligaments ensure the bone to bone connection both in guiding and in stabilizing the anatomical motion of the joints, keeping a physiological contact pressure with the articular surfaces (Benjamin & Ralphs 2000). Tendons transfer the force developed during muscle contraction to bones, allowing locomotion in the skeletal structure and enhancing joint stability (Benjamin *et al*. 2008). Since the mechanical function of ligaments and tendons is mainly related to bearing and to transmitting high tensile loads along their longitudinal axis, they have similar structural characteristics in order to confer a greater stiffness and resistance in the axial direction. The highly nonlinear and anisotropic strain–stress response of the tissues is due to the multi-level hierarchical organization of collagen type I network (Silver *et al*. 2003). Collagen type I represents the most abundant extracellular matrix (ECM) protein, and it has a well-defined axial orientation made of bundles embedded in a ground substance. The structural integrity and the viscoelastic characteristics for connective tissues result mainly from the interaction between collagenous proteins and non-collagenous proteins, e.g. proteoglycans (PGs; Cheng & Screen 2007). The microstructural arrangement of the ECM network determines the dynamical properties of force transmission inside the tissues, regulating the conversion of the strain energy and the performance of motor activities (Solomonow in press).

The interrelationships between the tensile characteristics of the tissue and the structural interactions of the ECM components at a microscopic level have attracted a broad range of interest in biomechanical and biological research. Classical models have addressed the passive structural behaviour of ligaments and tendons by combining a phenomenological approach (Kastelic *et al*. 1980; Butler *et al*. 1984; Danto & Woo 1993), with lower scale mechanisms of deformation (Diamant *et al*. 1972). Theoretical predictions have focused on quantifying relevant biomechanical parameters, such as the failure behaviour (Schechtman & Bader 2002), the hysteretic dissipation in cyclic tension (Maganaris 2002), the *in vivo* force measurement (Meyer *et al*. 2004) and the frequency spectrum of the strain-rate-dependent tensile properties (Pioletti & Rakotomanana 2000; De Vita & Slaughter 2006).

Molecular models have been used in computational biomechanics to confirm the ability of the discontinuous assembly of collagen type I and glycoproteins to transfer forces at the fibrillar level (Redaelli *et al*. 2003; Vesentini *et al*. 2005). Injury, in response to overloading, has been modelled as a rate-independent and irreversible process of structural damage occurring in collagen bundles and ground matrix (Gasser & Holzapfel 2002; Natali *et al*. 2004; Alastrue *et al*. 2007). Experimental studies have reported that the viscoelastic properties can be related to the mechanisms of strain transfer within the tissue, being dependent on the structural organization of the ECM components (Thornton *et al*. 1997, 2002; Screen 2007). The dissipative behaviour of PGs in the stress transfer between fibrils has been considered in microstructure-based finite viscoelastic models with the goal to reproduce the nonlinear softening in the mechanical response (Puxkandl *et al*. 2002). Creep, relaxation and hysteresis phenomena, finally, have been modelled from the observed microstructural deformation mechanisms at the fibrillar scale (Ciarletta *et al*. 2006, 2008).

The experimental investigation of tendons and ligaments has not defined a clear understanding of their microstructure, and an open debate concerns issues such as the fibril length and the fibrillar interactions (Scott *et al*. 2006; Vesentini *et al*. 2006). A multiscale modelling of the functional properties of the tissues in mechanobiology is needed to integrate the structural biology and the biochemistry into the biomechanical analysis of tissue behaviour (Wang 2006).

The aim of this work is to define a finite dissipative constitutive theory for ligaments and tendons based on the rate-dependent, strain-driven adaptive behaviour of the interfibrillar links between collagen type I fibrils. Outcomes from such a biomechanical approach can feed into experimental findings and may serve to support the interpretation of experimental data. In §2, we describe the biological characteristics of the ECM components in the fibrillar network, analysing the mechanical properties of collagen aggregates at different scales. The structural properties of PGs-mediated interfibrillar connections are discussed on the basis of the biochemical morphological evidence. The finite dissipative constitutive theory of ligaments and tendons is derived in §3 from thermodynamical considerations in a continuum mechanics framework. The viscoelastic response of the tissue at the macroscopic level reflects the thermally activated, strain-driven adaptive properties of fibrillar links by means of the theory of temporary networks. In §4, the proposed constitutive theory is applied to model the viscoelastic response of the tissues to cyclic uniaxial deformations. The mechanical behaviour under uniaxial tensile loading is analysed in §4.1, defining two energy-driven, strain-rate-dependent damage criteria for collagen fibrils and ground substance, respectively. The softening and the viscous dissipation during the tensile unloading are modelled in §4.2, using the concept of cooperative recruitment in transition state theory. The predicted stress–stretch responses of the proposed model are compared with the experimental behaviours on selected ligaments and tendons from referenced studies. In §5, we analyse the main results of the proposed dissipative model, evidencing the novelties and the limitations of the constitutive approach. Several applications are discussed in the wide research area ranging from medical practice to bioengineering science.

## 2. Microstructure and development of the mechanical properties in ligaments and tendons

### 2.1. Hierarchical organization of fibrillar collagen

The mechanical behaviour of dense connective tissues is primarily determined by the composition and the hierarchical organization of collagen molecules. Type I collagen is the basic and the most abundant ECM component in ligaments and tendons (Silver & Trelstad 1980). The monomer is a 300 nm long helix of three polypeptides, two α_{1} chains associated with one α_{2} chain, exceeding 1000 amino acid residues and obtained from the synthesized procollagen molecules through the cleavage of the terminal propeptides (Sun *et al*. 2002). The activation of the enzyme lysyl oxidase on the terminal telopeptides establishes intermolecular collagen cross links driving to a spontaneous linear aggregation in ordered supramolecular structures (Canty & Kadler 2002). The collagen molecules are staggered in the axial direction by approximatively 234 amino acid residues with respect to each other and form the fundamental structural unit called a fibril, having the characteristic diffraction pattern with a *D*-period of 67 nm (Hodge & Petruska 1963), with an overlap zone of 0.46*D* and a gap zone of 0.54*D* (Wess *et al*. 1998).

The collagen fibril grows by lateral aggregation of early subunits of few nanometres in diameter (Christiansen *et al*. 2000) and up to 100 μm in length (Birk *et al*. 1989), forming an ordered structure that has been reported to range from 20 to several hundred nanometres in diameter (Parkinson *et al*. 1997). The issue of fibril length in mature tissues is still controversial. Even if some authors suggest that fibrils may span the entire length of the tissue or may be considered functionally continuous, from a mechanical viewpoint a composite model of discontinuous fibril arrangement, up to several hundred micrometres in length, is better suited to explain the reported diffraction observations under tensile experiments (Vanderby & Provenzano 2003).

Few models have been presented to describe the spatial distribution of collagen molecules in fibrils (Prockop & Fertala 1998; Ottani *et al*. 2002). Even if a definitive scheme is far from being universally accepted (Koob & Summers 2002), a triclinic unit cell, where the molecules lie at a slight angle close to 5° to the fibril axis, has been identified from the crystallographic analysis (Ottani *et al*. 2001), thus suggesting a quasi-hexagonal lattice distribution (Orgel *et al*. 2001).

The hierarchical structure of collagen in ligaments and tendons results from the formation of a fine sheath of a loose connective tissue (Atkinson *et al*. 1999), the endotenon that provides the vascular, the lymphatic and the nervous supply inside the tissue. The endotenon enfolds bunches of parallel fibrils into fibres, in diameters ranging from 1 to 300 μm (Silver *et al*. 2003), and binds the fibres together in fascicles, as macroscopically visible in the section of the tissue (Yamamoto *et al*. 1999). Fibroblasts are sparsely aligned in rows along collagen fibres, having a characteristic spindle-shaped morphology and the function of synthesizing ECM proteins for healing and remodelling (Tohyama & Yasuda 2000). At the outer surface of the tissue, the endotenon continuously extends in order to form a loose connective layer of variable prominence, the epitenon, allowing the gliding function. Collagen type III is mainly found at the endotenon and the epitenon, with an expression ratio 1 : 9 compared with collagen type I, forming a smaller and less ordered fibrillar structure. The collagen organization in ligaments and tendons becomes fibrocartilagineous or fibrous in correspondence to the insertion sites to the bone, the entheses, having the function to provide structural continuity. The transmission of muscle-generated forces to tendons is ensured by the attachment of collagen fibrils to terminal myofibrils in the myotendineous junctions, whose stability is due to the highly infolded shape of the muscle cell membranes (Benjamin & Ralphs 2000).

### 2.2. Proteoglycans: biological function and interactions with fibrillar collagen

Proteoglycans (PGs) represent less than 1 per cent of the dry weight in tendons, being usually divided into two categories defined by the molecular size of the polypeptide core to which glycan chains are attached (Vogel & Heinegård 1985). Another classification separates two populations, PG I and PG II, based on different electrophoretic mobility and protease sensitivity. Ligaments and tendons are characterized by a predominant composition of small PGs of the PG II type, consisting of a core globular protein containing leucine-rich repeats (LRRs), and covalently bound to a maximum of two glycan chains (Vogel & Evanko 1987). Each LRR consists of 24 amino acids formed by leucine residues in fixed positions, and four species have been determined on similarities on the length and the sequence of the main stream.

Decorin and biglycan belong to a dominant subgroup in fibrillar matrix, showing 57 per cent protein sequence identity and containing either chondroitin sulphate (CS) or dermatan sulphate (DS; Svensson *et al*. 2001). The core proteins of these two PGs are of similar size, approximately 40 kDa, but they differ in the composition of the amino termini and, thus, in their functional roles (Järveläinen *et al*. 1991). Biglycan is biologically characterized by the presence of two glycosaminoglycan (GAG) chains and a greater molecular mass compared with decorin. It does not show any affinity with collagen type I, but it is localized close to specific cells implied in the regulation of their motility and adhesion processes. On the contrary, the synthesis of decorin by vascular cells is thought to be correlated to the gene expression of collagen type I (Nelimarkka *et al*. 1997). Other important features of decorin in tendons concern the binding capability with low-density lipoproteins and other collagen types (Pentikäinen *et al*. 1997) together with a regulation task in cell growth and migration (Svensson *et al*. 2001).

The core protein of decorin has been described as an arch-shaped protein formed by a central domain of 10 LRRs, having three attachment sites for N-linked oligosaccharides, whose three-dimensional structure consists of an inner concave surface covered by a curved β-sheet and an outer convex surface of α-helices (Weber *et al*. 1996). In this proposed geometry, the distance between the two arms of the arch is 6.5 nm, while the height and the overall thickness are 4.5 and 3 nm respectively, with two binding sites for fibrillar collagen located at LRRs 4–5 and composed of 40 amino acid residues (Svensson *et al*. 1995). While such a geometry would allow a collagen type I helix (diameter 1.5 nm) to be linked to the inner surface of the core protein (diameter 2.5 nm), a stable dimeric structure for the core protein of decorin in solution has been demonstrated, leading to a far less curved shape that suggests other ligand-binding scenarios for the fibrillar interaction (Scott *et al*. 2004).

Decorin links to fibrillar segments at the *d* band in the gap region during fibrillogenesis, functioning as a regulator of the ordered lateral growth and of the diameter and shape of the mature fibril (Pins *et al*. 1997). The GAG chains, having a length in the range of 60–160 nm, are highly negatively charged forming a typical bottlebrush-like structure with a relevant compressive strength. The arrangement of the GAG chains attached to the decorin core proteins provides the longitudinal stability of the interfibrillar structure in ligaments and tendons (Derwin *et al*. 2001). The rupture force of single bonds between two decorin core proteins was determined using laser tweezers as approximately 16.5 pN, supporting the hypothesis of a direct GAG chain interaction between decorins (Liu *et al*. 2005).

The GAG chain filaments of decorin have been visualized either orthogonal or parallel to the fibrillar surface in tendons (Raspanti *et al*. 1997). Orthogonal ridges encircle the entire fibrils in a regular array every D period; the length of these fragments is approximately the same as the centre-to-centre distance (approx. 66 nm) between fibrils (Scott 2001). Thus, the orthogonal CS/DS chains are supposed to interconnect adjacent fibrils in a regular network and to determine the ordered interfibrillar separation (Scott 1991). A skewness angle for these GAG chains has been observed increasing with the applied load, demonstrating the occurrence of a slippage during fibril-to-fibril interactions in loading (Liao & Vesely 2007). Peripheral axial filaments are sparsely detectable on the fibrillar surface, showing a length of approximately 1 D unit and a still not well-defined functional role (Scott 1988).

Minor PGs belonging to the same LRR family are fibromodulin and lumican. They show an interaction with collagen type I fibrils within the same binding site, which is different from that of decorin, and seem to be involved in the assembly of fibrillar collagen (Svensson *et al*. 2000). Another minor PG, the aggrecan, is found abundantly in fibrocartilagineous junctions for its ability to hold water and to resist compressive forces (Canty & Kadler 2002).

### 2.3. Mechanical properties of multiscale collagen aggregates

The mechanical properties of collagen aggregates show a wide variability from the molecular to the various hierarchical levels of organization in ligaments and tendons.

The collagen type I monomer is a non-rigid molecule assembled in consecutive regions of different flexibility: examinations by X-ray diffraction have evidenced linear elastic properties for the triple helix, with values of the elastic modulus ranging from 2.9 to 5.1 GPa (Sasaki & Odajima 1996*a*
; Sun *et al*. 2002). With an analogous study, a more compliant elastic behaviour for fibrils from bovine Achilles tendons has been reported, finding an elastic modulus of approximately 430 MPa and three possible deformation mechanisms for the staggered collagen molecules (Sasaki & Odajima 1996*b*
). Experimental data on collagen fibres show a non-Hookean elastic behaviour characterized by a stress–strain curve with an initial toe region and a subsequent linear region, with elastic modulus ranging from 269 to 484 MPa (Gentleman *et al*. 2003).

Recent studies have demonstrated that the geometry of the fibrillar assembly strongly affects the mechanical properties at a macroscopic level (Christiansen *et al*. 2000). Both the ultimate tensile strength and the elastic modulus are linearly dependent on the fibril length, while a greater fibrillar diameter increases the resistance to deformations only at low strain levels (Silver *et al*. 2003). In addition, the presence of decorin increases the tensile properties of collagen fibres allowing a more efficient slippage of the fibrillar units (Pins *et al*. 1997).

For tendon fascicles, the tensile modulus increases nonlinearly for decreasing cross-sectional areas of the samples, with an average value of 200 MPa (Atkinson *et al*. 1999).

The fibrillar structure has been investigated in the submicrometre range using synchrotron X-ray scattering during tensile loading of the macroscopic tissue (Fratzl *et al*. 1998). The mechanical stress–strain curve begins with a characteristic toe region, with a visible incremental straightening of crimped collagen fibrils, and a progressive disappearing of the vertical crimp banding (Hansen *et al*. 2000, 2002). For strains beyond 3 per cent, typical increase in stiffness of the heel region is associated with the occurrence of molecular kinks by thermal activation. At a higher level of strain, approximately 5 per cent, entropic elasticity is progressively replaced by the linear elastic stretching of the collagen type I triple helices, corresponding to an increase in the D band at the light microscope (Mosler *et al*. 1985).

The functional adaptation of ligaments and tendons results in a reduction of shear modulus, tensile strength and ultimate strain in consequence of stress shielding (Yamamoto *et al*. 1999). The biomechanical remodelling is ensured by the role of fibroblasts in the ECM, whose presence creates a local stress enhancement at the interface with fibrillar collagen and prevents a detrimental effect on the entire tissue (Tohyama & Yasuda 2000).

## 3. Dissipative constitutive equation for the anisotropic response of ligaments and tendons with rearranging ECM network

Ligaments and tendons can be modelled as multi-level composite tissues, made of a directionally oriented network of discontinuous collagen fibrils immersed in a hydrated matrix. The high water content that characterizes the tissues allows one to consider them as incompressible continuum bodies, whose time-dependent mechanical behaviour is primarily due to structural modifications occurring in the fibrillar network (Cheng & Screen 2007).

Let *Ω*
_{0} be the fixed reference configuration of the soft collagenous tissue in the unloaded state. The description of the deformation can be defined by an objective mapping *Χ*: *Ω*
_{0}→
^{3} that transforms the material point **
X
**∈

*Ω*

_{0}to a position

**=**

*x**Χ*(

**,**

*X**t*) in the deformed configuration

*Ω*. A time-dependent strain-energy function

*Ψ*can be introduced as a function of the relative deformation tensor (

*t*,

*τ*)=Grad(

*t*,

*τ*)

*Χ*(

**,**

*X**t*), where Grad(

*t*,

*τ*) indicates the gradient in material coordinates between the final configuration at time

*t*and the reference configuration at time

*τ*. Such a viscoelastic strain-energy function

*Ψ*must obey the principle of material frame indifference, being invariant for any superimposed rigid-body motion 3.1 being and ; where indicates the space of linear transformations in .

Strain measurements, obtained using local markers on a confocal microscope, have demonstrated that the effective stretch of the collagen fibrils inside ligaments and tendons is much less than that of the macroscopic tissue (Screen *et al*. 2004). This local strain field within the tissue is dominated by a sliding behaviour between fibrils, which provides the major extension mechanisms for the overall deformation. The PGs-mediated interactions between fibrils are postulated to be the viscous mechanical structures that transmit such a considerable deformation outside the fibrillar collagen. Recent biomechanical models have demonstrated through a multiscale approach that GAG interfibrillar bridges associated with decorin may provide structural integrity to the tissue, being able to transfer forces between discontinuous fibrils (Puxkandl *et al*. 2002; Redaelli *et al*. 2003). Even if the binding configuration between the dimeric structure of decorin and the fibrillar collagen is still unclear, a computational analysis accounting for a monomeric decorin structure has indicated that such a fibril-to-fibril assembly implies GAG chain rupture in response to overload (Vesentini *et al*. 2005). Moreover, experimental observations of the tertiary structure of GAG chains have evidenced a stress-driven reversible rupture and cooperative slippage of GAG-mediated junctions, with the occurrence of irreversible strains if overlaps between chains decrease beyond a critical level (Scott 2003).

In order to account for the observed microstructural interactions, continuous breakage and reformation processes of the adaptive GAG-mediated junctions need to be considered in order to derive a dissipative constitutive model of ligaments and tendons (Ciarletta *et al*. 2006).

Let *N*
_{m}(*t*, *τ*), *N*
_{fj}(*t*, *τ*) be the number per mass unit of the isotropic matrix and of the generic collagen fibrils with preferred *j* direction, respectively, considered rearranged before the intermediate time *τ* and still active at the actual time *t*, being 0≤*τ*≤*t*.

Considering the transient network theory, the relative rates of reformation *γ*
_{m}(*τ*), *γ*
_{fj}(*τ*) of the rearranging matrix and fibrillar units, respectively, are modelled to occur at random times and defined as (Drozdov 1998, 2000)
3.2
Similarly, the rates of breakage of the active regions of matrix and collagen fibres *Γ*
_{i}(*t*,0), *Λ*
_{i}(*t*, *τ*) can be expressed as
3.3
3.4
The time-dependent strain-energy function for the soft collagenous tissue can be decoupled as the sum of an isotropic term *Ψ*
_{m}(
, *t*) and of the anisotropic terms *Ψ*
_{fj}(
, *t*), representing the contribution of the hydrated matrix and of the collagen fibrils in the *j* direction, which are active at time *t* (Gasser *et al*. 2006):
3.5
where *p* is the Lagrange multiplier introduced to restore the incompressibility condition.

The initial number densities *N*
_{m}(0, 0), *N*
_{fj}(0, 0) for the matrix and the collagen fibrils are assumed to be independent of the reformation process for a preconditioned tissue. Under these assumptions, the strain-energy functions *Ψ*
_{m}(
, *t*) and *Ψ*
_{fj}(
, *t*) can be defined as
3.6
where *ψ*
_{i} correspond to hyperelastic strain-energy functions and
(*t*, *τ*) is the deformation tensor for a portion rearranged at time *τ* and active at time *t*.

Equations (3.5) and (3.6) define a time-dependent strain-energy function for a soft collagenous tissue considered as a composite material with adaptive regions that fail and rearrange at random times, accounting both for recoverable and irreversible GAG chain detachments.

The terms *ψ*
_{i} in equation (3.6) define the hyperelastic energy associated to the tissue response, considering that the external loading is applied at sufficiently high strain rate so that the energy dissipation can be neglected. The isotropic hyperelastic contribution of the ground matrix in equation (3.6) can be expressed by a Mooney–Rivlin constitutive equation
3.7
where *I*
_{k} represent the *k*th invariants (*k*=1, 2) of the right Cauchy deformation tensor
, while *c*
_{1} and *c*
_{2} are material constants of the tissue whose values depend on the content of ground substance in the tissue.

For what concerns the anisotropic behaviour of the tissue, let us consider a single preferred material direction of collagen fibre reinforcement, described by the unit vector **
j
**. In the case of more complex collagen structures, such as planar collagenous tissues (oesophagus, intestine, blood vessels), an extension including more preferred directions is straightforward. Recalling equation (3.1), the principle of material frame indifference for the directional reinforcement requires the definition of the anisotropic hyperelastic strain energy as a function of the two independent pseudo-invariants that result from the combination of
, and the structural tensor

**⊗**

*j***: 3.8 The pseudo-invariant**

*j**I*

_{5j}expresses the kinematical behaviour of fibrils in correspondence with shear strain states. As the cross-sectional radius of the collagen fibrils is very small compared with their length, the influence of

*I*

_{5j}can be included into a functional dependence of the anisotropic strain-energy function on

*I*

_{4j}as follows: 3.9 where

*k*

_{1j}is a constant stress-like material parameter related to the concentration of collagen fibrils, and

*k*

_{2j}is a dimensionless material constant related to the crimping level of their initial conformation. Even if other expressions for the anisotropic strain-energy function may be considered (Holzapfel

*et al*. 2000; Gasser

*et al*. 2006), the choice of equation (3.9) is done to allow a direct comparison of the numerical results of this model with those from previous studies on ligaments and tendons (Natali

*et al*. 2004, 2005).

In order to derive a constitutive equation for the soft tissue, with the strain-energy function as defined by equations (3.5)–(3.7) and (3.9), let us consider the expression of its time derivative:
3.10
The reformed matter is considered in the stress-free configuration at the instant of its creation, so that
(*t*, *t*)=0. The term *Φ*(*t*) in equation (3.10) is defined as
3.11
Recalling that the rearranging processes of the fibrillar units and of the hydrated matrix are assumed to occur at random times, the differentiation of the tensor
(*t*, *τ*) leads to the following tensorial relationship:
3.12
where
is the rate of strain tensor at time *t*, being
(*t*, 0)=
(*t*, *τ*).

Considering a preconditioned tissue, from equation (3.3) the number *N*
_{i}(*t*, 0) of the active fibrils at time *t* can be expressed as function of its initial number density of active matter as
3.13
Under the same assumption, the number of the rearranged regions per unit mass is given, from equation (3.4), by the following relationship:
3.14
The chain rule for differentiation of the hyperelastic energies leads to the following:
3.15
3.16
being
3.17
3.18
Let us recall the following tensorial properties:
3.19
and the pull-back operation of the rate of strain tensor
in terms of the time derivative of the Cauchy–Green strain tensor
, being
. The expression in equation (3.10) of the time derivative of the strain-energy function can be rewritten as
3.20
being
3.21
3.22
where **
b
**(

*t*,

*τ*) indicates the left Cauchy strain tensor for a portion rearranged at time

*τ*and still active at time

*t*.

In order to derive a constitutive equation, let us consider the Clausius–Planck form of the second law of thermodynamics for an isothermal process (Holzapfel 2001):
3.23
where
is the second Piola–Kirchhoff stress tensor (the work conjugate of
), and *Π*
_{int} represents the internal dissipation or the local production of entropy.

After simple manipulations, it can be demonstrated that the term *Φ*(*t*) in equation (3.11) fulfils the following inequality:
3.24
if *Λ*
_{i}(*t*, *τ*), *γ*
_{i}(*τ*) are positive definite functions for 0≤*τ*≤*t*, being *i*=(*m*, *fj*).

Considering the inequality in equation (3.24), a dissipative constitutive equation for isothermal loading can be derived equating to zero the term between square brackets in equation (3.23):
3.25
In the constitutive theory expressed by equation (3.25), the local production of entropy is entirely due to the term *Φ*(*t*), being related to the kinetics of the rearranging processes as defined in equation (3.24).

The constitutive relationship in equation (3.25) is fully determined by the isotropic hyperelastic energy of the matrix, in equation (3.7), by the anisotropic hyperelastic contribution of the collagen fibrils, in equation (3.9), and by the thermally activated processes that control the kinetics of breakage and rearrangement inside the tissue, determining the expressions of *Λ*
_{i}(*t*, *τ*), *γ*
_{i}(*τ*), at time 0≤*τ*≤*t*.

## 4. Analysis of damage, softening and non-recoverable strain in the hysteretic spectrum of ligaments and tendons

The aim of this section is to apply the proposed constitutive theory of adaptive GAG-mediated junctions between collagen fibrils to study the functional adaptability of ligaments and tendons. In particular, the study will focus on the following characteristics of the mechanical response under a cyclic external strain:

— the occurrence of a damage and its rate-dependent laws of evolution during the loading phase until rupture,

— the effect of the stress-softening phenomenon, the so-called Mullins effect, on the viscoelastic energy dissipation,

— the presence of a residual strain after the application of a cyclic external deformation, and

— the influence of the stress memory of the tissue in determining its hysteretic behaviour.

Recalling the applicability of the proposed model to a generic multiaxial state of stress, the mechanical behaviour of ligaments and tendons is considered under uniaxial tensile tests in the following. Owing to the ECM structural organization of the tissues and their physiological functions, this approach allows an immediate physical interpretation of the results and an easy correlation with experimental data.

In order to test the predictive ability of the proposed dissipative constitutive theory, the analytical model of adaptive GAG-mediated junctions is compared with several experimental stress–stretch curves from the biomechanical literature. The procedure of curve fitting with the experimental data is performed using the Levenberg–Marquardt algorithm in order to determine the mechanical parameters. The numerical integrations are calculated using recursive Lobatto quadrature. An optimization method is implemented to solve in the least-square sense the nonlinear fitting of multiple experimental curves. Extended details about the experimental methods can be found in the referenced papers.

### 4.1. Strain-driven and strain-rate-dependent evolution equations of damage in the loading phase

Classical hyperelastic models, as the ones expressed by the constitutive relationships in equations (3.7) and (3.9), aim at describing the mechanical response of the macroscopic tissue during normal physiological functions, having a biomechanical validity that is restricted to a deformation range where no damage has been induced to the material.

In general, the evolution of a damage occurring inside a soft tissue is a discontinuous phenomenon possessing a well-defined activation energy. In the following, we will consider the damage for the matrix and the collagen fibrils, separately, applying the Lemaitre–Chaboche formalism to describe the failure of the components of the adaptive ECM network (Lemaitre & Chaboche 1985).

Let us perform a uniaxial tensile test, with a generic strain rate *v*, over the principal direction having unit vector **
e
**

_{1}=

**. The symmetry of the problem allows us to consider a transversely isotropic behaviour of the tissue, being**

*j**λ*(

*t*)=(1+

*vt*) the time evolution of the principal stretch. Considering the general incompressibility condition, the principal stretches

*λ*

_{i}are linked by the simple relationship Considering a traction uniaxial tensile test starting at

*t*=0, the value of the pseudo-invariant corresponding to the maximum collagen fibre stretch, measured along the principal direction

**and up to the time**

*j**t*, is expressed by 4.1 A damage surface in the strain space for the collagen fibres can be defined by the following relationship (Rodriguez

*et al*. 2008): 4.2 In order to define the evolution process of the damage, the normal

*n*_{fj}to the damage surface can be written as 4.3 Since

*δ*is an arbitrary admissible variation of the right Cauchy tensor, damage evolution occurs if and only if the following condition holds (Balzani

*et al*. 2006): 4.4 Considering the anisotropic damage phenomenon in ligaments and tendons as a breakage process of GAG-mediated interfibrillar links that cannot be recovered, the strain-energy function in equation (3.9) is related uniquely to the number

*N*

_{fj}(

*t*,0) of collagen fibrils per mass unit initially existing at time

*t*=0, and still active at time

*t*.

The anisotropic damage during the loading phase can be described, from equation (3.22), by the following constitutive equation:
4.5
From equations (3.3) and (4.4), the evolution of the anisotropic damage can be expressed as
4.6
Following a similar approach, a damage surface for the isotropic matrix can be identified as follows (Chagnon *et al*. 2004):
4.7
*N*
_{m}(*t*, 0) being the number density of hydrated matrix initially existing at time *t*=0 and still active at time *t*, the constitutive equations of the isotropic damage evolution are the following:
4.8
4.9
The Eyring theory of thermally activated processes can be used to define the rates of breakage of the active regions of matrix and collagen fibrils (Eyring 1936), *Γ*
_{i}(*t*, 0) being defined as
4.10
where
are the constant rates of breakage at the occurrence of damage, and
is the excess energy that must be accounted to include damage evolution. According to the damage evolution equations expressed by equations (4.6) and (4.9), such excess energies
can be defined as
4.11
Equation (4.11) states that the excess energies for the matrix and the collagen units grow slower than their respective hyperelastic energy, with a dependence expressed by the Wagner phenomenological relationship, *κ*
_{i} being constant material parameters (Wagner & Schaeffer 1992).

Considering the constitutive relationships in equations (3.25), (4.5) and (4.8), and deriving the hydrostatic pressure *p* from the boundary conditions *S*
_{22}(*t*)=*S*
_{33}(*t*)=0, the principal component *S*
_{11}(*t*) of the second Piola–Kirchhoff tensor in the direction of collagen reinforcement can be expressed as
4.12
the terms *η*
_{m}(*t*), *η*
_{fj}(*t*) being, in accordance with the damage evolution laws in equations (4.6) and (4.9), expressed by the following:
4.13
where *λ*
^{pr} represents the maximum stretch level reached during the preconditioning of the tissue, being
.

Equations (4.12) and (4.13) are the constitutive relationships defining the time-dependent uniaxial tensile response of ligaments and tendons accounting for the energy-driven damage of fibrillar collagen and of hydrated matrix.

The predicted uniaxial loading curve for ligaments and tendons has a typical failure behaviour characterized by two stress peaks with a strain-rate-independent ultimate stretch (Wren *et al*. 1998; Liao & Belkoff 1999), corresponding to progressive failure of collagen fibrils and of the isotropic ECM, as shown in figure 1. The rupture of the ECM corresponds to a lower stress peak with a higher ultimate stretch, occurring when almost all fibrils have failed.

The predictions of the constitutive model of strain-driven damage evolution in equation (4.12) have been compared with two sets of experimental data on extensor digitorum longus tendons (EDLT) reported in Schechtman & Bader (1997), and the results are shown in figure 2. The curves in figure 2 show the predictive ability of the theoretical model but they do not allow a direct comparison of the damage behaviours: it is well known that the failure response of the tissue strongly depends on the relationship between the width and the length of the sample (Zienkiewicz & Taylor 1991).

Considering a uniaxial loading with a high strain rate, the damage evolution laws can be considered as the dominant viscous process inside the tissue. The viscoelastic response is characterized by the strain-rate-dependent damage evolution laws in equation (4.13). The loading behaviour in equation (4.12) has been compared with the strain-rate-dependent stress data of patellar tendon (PT) and of anterior cruciate ligament (ACL), from Pioletti *et al*. (1998) and Pioletti & Rakotomanana (2000), as depicted in figure 3.

### 4.2. Cooperative relaxation and rearrangement of the adaptive ECM network during the unloading phase

Let us consider the modelling of the adaptive mechanical response of ligaments and tendons during the unloading phase of the uniaxial tensile test at strain rate *v*. If *t*
^{*} is the time corresponding to the maximum principal stretch, *λ*
^{*}=*λ*(*t*
^{*}), then the deformation tensors can be written as
4.14
where *λ*(*s*)=*λ*
^{*}−*v*.*s*, within the range 0≤*s*≤(*λ*
^{*}−1)/*v*.

The viscoelastic tensile properties of ligaments and tendons are characterized by varying behaviour of different ECM portions ranging from a crystalline conformation to a state resembling a liquid crystal. The existence of a shrinkage transition has been demonstrated by measuring the potential energy barrier that should be overcome by a tissue region to pass into to the relaxed state (Cohen *et al*. 1976). This potential energy has been incorporated into a viscoelastic constitutive equation at finite deformation to model the relaxation and creep spectra in flexor digital tendons (FDT; Ciarletta *et al*. 2006). More recently, the adaptive fibrillar recruitment into fluid-like states has been related to the softening effect in the tensile behaviour of tendons during the unloading phase, characterizing their experimental hysteretic spectrum at physiological deformations (Ciarletta *et al*. 2008).

These experimental observations lead one to consider the cooperative relaxation of adaptive portions of collagen fibres and hydrated matrix as a thermally activated process depending on the hyperelastic strain-energy function at a generic time *t*. From the kinetic relationships in equations (3.3) and (3.4), we assume that
4.15
4.16
As described by evolution laws in equations (4.6) and (4.9), a damage inside the tissue cannot be generated during the unloading phase. As a consequence, the total number per unit mass of existing rearranging collagen and matrix portions is constant at each *t*
^{*}≤*t*(λ^{*}−1)/*v*, imposing the following relationship:
4.17
*N*
_{i}(*t*, *t*)=*N*
_{i}(*t*
^{*},0)=*η*
_{i}(*t*
^{*}).*N*
_{i}(0, 0) being the constant number density of rearranging portions.

In order to calculate the mechanical response of the collagen reinforcement during the unloading case, we should consider that the rearranged fibres do not actively contribute, being *I*
_{4m}(*t*, *τ*)=[*λ*(*t*)/*λ*(*τ*)]^{2}≤1 for each *τ*≤*t*. The anisotropic tensile response
can be finally expressed as
4.18
being
4.19
The isotropic mechanical response
during the unloading phase is composed by the sum of tensile contribution due the material already existing at time *t*
^{*}, and the compressive contribution due to the portions of the hydrated matrix that have been reformed at a time *τ*>*t*
^{*}. From the constitutive relationship in equation (3.21), the analytical expression of
is the following:
4.20
Recovering the incompressibility constraint *p* from the boundary conditions
, the principal component
of the second Piola–Kirchhoff tensor in the unloading phase can be expressed as
4.21
Equation (4.21) describes the unloading tensile response of the tissue including a softening effect resulting from the cooperative relaxation of the ECM components, and a residual strain due to the strain-activated reformation process inside the hydrated matrix. In order to determine the time evolution of the term under the integral sign in equation (4.21), the differentiation of equation (4.17) leads to the following differential equation:
4.22
which is valid for a generic function *f*[*λ*(*τ*)].

The hysteretic response of ligaments and tendons at physiological strain levels is characterized by a quasi-flat hysteretic spectrum at high stretch limits, with a dispersed energy per loading–unloading cycle decreasing with an increasing strain rate. The dissipative constitutive response in the unloading phase, determined by equation (4.21), is compared with the experimental data on FDT, from Ciarletta *et al*. (2008). The theoretical prediction of the proposed model is shown in figure 4 for different values of the maximum stretch level in the load cycle.

Ligaments and tendons act as optimized mechanical buffers in high-frequency actions at physiological deformations, keeping a memory of the structural damage accumulated in their loading history. The predictions of the hysteretic response of the proposed dissipative model are compared with the experimental cyclic response of the periodontal ligament (PDL), from Natali *et al*. (2004). The loading and unloading responses of the dissipative model, from equations (4.12) and (4.21), are shown in figure 5 at different strain-rate values of the uniaxial traction cycle.

## 5. Discussion and concluding remarks

The authors have proposed in this study a nonlinear dissipative constitutive theory for ligaments and tendons based on thermodynamic considerations on the multiscale structure of their temporary ECM network. The viscous biomechanical behaviour of the anisotropic network of discontinuous collagen type I fibrils has been modelled from the adaptive properties of the PGs-mediated interfibrillar junctions, considered embedded in an isotropic hydrated matrix. The novelty of this research is to include both reversible (slip-links) and irreversible (rupture mechanism) detachment of the GAG junctions between collagen type I fibrils with the aim to describe the thermally activated, strain-driven behaviour of the tertiary structure of the GAG chains of decorin in fibril-to-fibril connections. The purpose of this approach is to consider the effects in tissue viscoelasticity due to the unfolding mechanism preceding fibrillar elongation, together with strain-driven kinetic slippage of entire fibrils and the strain-rate-dependent damage characteristics of the interfibrillar bridges. The viscoelastic strain energy of the temporary ECM network is defined as the sum of the contributions of the isotropic ground substance and of the anisotropic fibrillar reinforcement, defining the number density of the active matter through the kinetic behaviour of thermally activated breakage and reformation processes. The constitutive theory is derived from the Clausius–Planck form of the second law of thermodynamics for an isothermal process, ensuring the positive definiteness of energy dissipation for the proposed model. The isotropic and the anisotropic damage of the tissue are considered using the theoretical framework of continuum damage mechanics in order to describe injury as the irreversible failure of the biological constituents of the temporary ECM network with different activation energies. The dissipative description of the damage inside ligaments and tendons is derived from the Eyring theory applied to the irreversible breakage process of active ECM matter. The transition state theory considers a strain-driven activation energy of the failure process resulting in two strain-rate-dependent damage evolution laws for fibrillar units and ground matrix, respectively. Typical tensile responses to uniaxial traction tests of ligaments and tendons including strain-dependent damage are depicted in figure 1, from the theoretical predictions of equations (4.12) and (4.13). Such theoretical predictions of the proposed model are consistent both with the typical tensile responses in the experimental post-failure behaviour of the tissues (Liao & Belkoff 1999) and with the strain-rate-dependent peak stress corresponding to an invariable ultimate stretch limit (Wren *et al*. 1998). The influence of preconditioning on the tensile response of the material is described as a spatial memory of the active fibrillar units resulting from the load history of the soft tissue (Zhu *et al*. 1996). The theoretical predictions of the model are compared with the experimental data of the tensile response under uniaxial loading of ligament and tendon samples in the biomechanical literature. The theoretical and experimental curves are shown in figures 2 and 3, while the fitted microstructural material parameters are collected in table 1, demonstrating the ability of the proposed model to describe accurately the viscoelastic damage behaviour of the tissues under the loading phase.

In order to model the hysteretic response of ligaments and tendons under a cyclic uniaxial strain, the softening of the macroscopic tissue in the unloading phase is related to the adaptive microstructural recruitment of the ECM constituents (Ciarletta *et al*. 2006). Reversible interactions of GAG chains are modelled as a strain-driven slippage mechanism causing a progressive detachment of fibrillar collagen, followed by a GAG bridging process from the resting state of the reformed interfibrillar connections (Scott 2003). A similar cooperative process is assumed for the temporary structure of the hydrated ground matrix, reflecting the energy-activated transition from a crystalline conformation to a liquid-like state observed experimentally (Prockop & Fertala 1998).

The time-dependent response of ligaments and tendons under a cyclic uniaxial tensile test is modelled by the constitutive relationships in equations (4.12) and (4.21) for loading and unloading, respectively. The proposed dissipative theory defines a fully nonlinear hysteretic behaviour from the microscopic deformation mechanisms of the temporary ECM network. Damage evolution, stress memory, softening phenomena and non-recoverable strain are modelled in equations (4.13) and (4.19) through the definition of material parameters representing the adaptive mechanical behaviour of the ECM components at the fibrillar level. The ability of the theoretical predictions to reproduce the experimental data of the complex time-dependent finite response of the tissues is shown in figures 4 and 5, with the fitted material parameters collected in table 2. The variability in the range of values for the mechanical parameters used in the theoretical predictions, as collected in tables 1 and 2, is due to the morphological and compositional differences between ligaments and tendons, and to the structural specificity of the tissues related to their mechanical function (Amiel *et al*. 1984; Benjamin & Ralphs 2000). In general, tendons have a higher content of collagen type I than ligaments, determining a variability in the parameters (*c*
_{1}, *c*
_{2}) and *κ*
_{1j}, with a variable undulation degree in the crimp pattern that influences the value of *κ*
_{2j}. Moreover, ligaments have a larger number of small fibrils and a higher density of PGs compared with tendons, resulting in different values of the mechanical parameters determining the damage
and the softening
characteristics of the tissue.

The proposed constitutive model accounts for the microstructural mechanisms of deformations of the complex PGs–collagen type I, investigating the biomechanical role of the ECM components in determining the dynamic tensile properties of ligaments and tendons (Puxkandl *et al*. 2002; Lucas *et al*. 2008). The results of the model are consistent with the experimental results obtained with transgenic mice with altered ECM: the proposed analysis predicts, in fact, that the GAG content is the dominant structural factor in determining the viscoelastic properties of the macroscopic tissue (Elliot *et al*. 2003). Moreover, the collagen type I content is found to have a direct influence on the overall stiffness, while a variation in the GAG content is observed to cause a dominant effect in determining the mechanical properties (Robinson *et al*. 2004). The softening effect due to the fibrillar sliding in the temporary network is confirmed by the experimental analysis of the rupture properties of tendon samples incubated in enzymatic solutions designed to disrupt or swell the non-collagenous matrix (Screen *et al*. 2006). In addition, experiments on decorin-deficient mice have demonstrated the key regulatory function of this PG to coordinate the process of fibrillogenesis and tissue development (Zhang *et al*. 2006).

A limitation of the proposed constitutive model is represented by the hypothesis of homogeneity of the fibrillar structure at the interface that could possibly mimic an altered adhesive behaviour of the tissue (Tang *et al*. 2005). The description of the microstructural ECM deformation should include the presence of sparsely oriented mineralized fibrils in entheses that have been reported to influence the failure behaviour for inhibiting collagen sliding (Siegmund *et al*. 2008). Future work should investigate the dissipative effects in a more complex morphology in the ECM network, accounting for the presence of helically oriented fibrils and GAG chains non-orthogonal to the collagen fibrils (Henninger *et al*. 2007). A non-isothermal constitutive theory may be derived in order to model the detrimental effect in the ECM structure induced by the variation of PH conditions and temperature (Usha & Ramasami 1999).

In conclusion, knowledge of the correlation of the viscoelastic behaviour with changes in the microstructure of ligaments and tendons is fundamental to optimize the design of functional engineered tissues for proper replacement or repair purposes (Ahearne *et al*. 2005). Culture conditions and mechanical stimuli applied to collagen matrices have been demonstrated to play a fundamental role in the synthesis and the assembly of the ECM components, inducing a strain-dependent expression of decorin in the cellularized matrix (Ferdous *et al*. 2008). The balance between the stiffness and the viscoelastic properties of new scaffolds may be obtained through an interplay of mechanical stimulation, growth factors and hormones for an optimal control of the tissue engineering process (Hoffmann & Gross 2007). Experimental results indicate that collagen synthesis is controlled by a complex mechano-sensitive process with a relevant temporal dependence (Maeda *et al*. 2007), suggesting that a microstructural biomechanical approach may be important for the understanding of the mechanotransduction paths in tissue remodelling and healing. Further applications of the proposed constitutive model can be envisaged, finally, in bioengineering science and medical practice, bridging together robotic technology and computational simulation for the design of implantable devices in prosthetics (Cavallaro *et al*. 2005) and the clinical treatment of injury (Woo *et al*. 2006).

## Footnotes

- Received November 14, 2008.
- Accepted November 25, 2008.

- © 2008 The Royal Society