## Abstract

Soft tissues exhibit complex viscoelastic behavior, including strain-rate dependence, hysteresis, and strain-dependent relaxation. In this paper, a model for soft tissue viscoelasticity is developed that captures all of these features and is based upon collagen recruitment, whereby fibrils contribute to tissue stiffness only when taut. We build upon existing recruitment models by additionally accounting for fibril creep and by explicitly modeling the contribution of the matrix to the overall tissue viscoelasticity. The fibrils and matrix are modeled as linear viscoelastic and each fibril has an associated critical strain (corresponding to its length) at which it becomes taut. The model is used to fit relaxation tests on three rat tail tendon fascicles and predict their response to cyclic loading. It is shown that all of these mechanical tests can be reproduced accurately with a single set of constitutive parameters, the only difference between each fascicle being the distribution of their fibril crimp lengths. By accounting for fibril creep, we are able to predict how the fibril length distribution of a fascicle changes over time under a given deformation. Furthermore, the phenomenon of strain-dependent relaxation is explained as arising from the competition between the fibril and matrix relaxation functions.

## 1 Introduction

Soft tissues that connect, support, or surround bones and organs are abundant. Such tissues are generally inhomogeneous and have multiple phases, which provide a challenge to those attempting to understand their mechanical behavior. The macroscopic response of a tissue to a force depends on the constitutive behavior of each of its individual phases, as well as the manner in which these phases interact.

Connective tissues are commonly fibrous, being composed of a compliant material within which a number of stiffer fibers reside. Tendons, for example, are composed of collagen fibrils embedded in a matrix that consists largely of proteoglycans and elastin [1]. The specific volume fraction, distribution, and arrangement of the fibers depend on the tissue type and often vary throughout the tissue.

In this paper, we develop a microstructural model of the macroscopic viscoelastic behavior of tendon or ligament fascicles. Although tendons and ligaments have different functions (tendons connect bone to muscle, whereas ligaments connect bone to bone), they are composed of the same fundamental components. Their main subunit is the fascicle, which consists of crimped collagen fibrils (see Fig. 1) that have a diameter of 50–500 nm [2]. Fibrils are in turn composed of microfibrils that are made from a regular arrangement of collagen molecules, but we shall not consider the mechanics of any lengthscales below the fibril level in this paper. The constitutive behavior of individual collagen fibrils has been reported to be linearly viscoelastic [35].

The goal of understanding the constitutive behavior of soft tissues, and more specifically here, tendons, and ligaments, is of fundamental importance. One area that could benefit from such an understanding is tendon reconstruction surgery. A replaced or repaired tendon must possess similar mechanical properties to the original tissue in order to carry out its role effectively [6]. Therefore, to design effective prosthetic tendons, we must have knowledge of the relationship between tissue microstructure and mechanical behavior.

Tendon and ligament fascicles in general exhibit strong nonlinear viscoelasticity with time-dependent and load-history dependent mechanical behavior even under small strains. The fact that the constitutive response of tendon and ligament is nonlinear and viscoelastic was established some time ago. For example, Rigby et al. measured the properties of rat tail tendons [7]. By viewing fascicles with polarized light, they attributed the nonlinearity to the successive straightening of crimped collagen fibrils.

A simple way to approximate the constitutive behavior of soft tissue is to ignore viscous effects and consider a purely elastic model. Many elastic models have been proposed that can be categorized according to their phenomenological or microstructural basis. The former approach often utilizes nonlinear elasticity to model soft tissues; Fung [8], for example, modeled rabbit mesentery constitutive behavior using an exponential function, and more recently, Holzapfel et al. [9] developed a strain energy function for the modeling of arteries that similarly utilizes an exponential function.

Microstructural models are appealing in the sense that the actual mechanics of each subunit can be modeled and upscaled; however, it is often unclear exactly how such microstructural elements behave and interact as their properties need to be determined experimentally, which can be a huge challenge in itself. A popular elastic model is that of Kastelic et al. [10], which introduced a mechanism to include the notion of sequential straightening and loading (SSL) to describe the uncrimping of fibrils referred to above. Prior to this, Lanir [11] had proposed a model which had some similarities, but which dealt with the uncrimping in a different way. He assumed that the crimp is induced and sustained by filaments that are attached to the fibril at random intervals; however, although interesting, this approach has not proved to be as popular. Ault and Hoffman [12,13] adapted micromechanical models originally devised for fiber reinforced composites to tendon, and Shearer [14] extended the SSL model to derive a new strain energy function that is valid for general deformations within a nonlinear elastic framework, and which is still based on tendon and ligament microstructure. This model was shown to fit experimental data accurately.

Some elastic models assume the fibrils are helically arranged: Beskos and Jenkins [15] considered tendons whose fibrils were assumed to be helical and inextensible, which led to their model having infinite stiffness at full extension. Freed and Doehring [16] and Grytz and Meschke [17] extended that work by modeling the fibrils as helical springs and Shearer [18] considered fibrils that are helical and crimped. Other authors, such as Hurschler et al. [19], assumed a statistical distribution of fibril alignments in order to determine their average stress–strain properties in a given direction.

A number of models of soft tissue viscoelasticity have also been proposed. Linear viscoelastic models were used early on [20] and have the advantage of being easy to implement mathematically. Viidik [21] later used an assembly of springs, dashpots, and slack elements to model the load-deformation curve of tendon. In reality, however, tendon viscoelastic behavior is nonlinear. Experiments show that rates of relaxation and creep are dependent on the strain or stress level that is being imposed [22,23]. The latter rules out the possibility of employing quasi-linear viscoelasticity (QLV) (which is based on the work of Fung [24] and has recently been reinterpreted by De Pascalis et al. [25], extended to the case of transverse isotropy by Balbi et al. [26] and employed for modeling viscoelastic inflation problems by De Pascalis et al. [27]) with a single scalar relaxation function. QLV assumes that the viscous relaxation rate is independent of the instantaneous local strain and is a special case of the more general constitutive model developed by Pipkin and Rogers [28]. One of the earliest examples of using QLV to model ligaments was by Woo et al. [29], who showed that QLV can achieve good agreement with individual experiments, but they did not consider strain-dependent relaxation.

Although QLV with a single scalar relaxation function gives different stress–strain curves for different strain rates, the fact that the relaxation function is independent of deformation means that an imposed step-function in strain leads to the same rate of relaxation regardless of the magnitude of the imposed strain. Such models have been fitted successfully to a number of datasets, although frequently these only incorporate data from experiments at a single strain level or strain rate [30]. On the other hand, DeFrate and Li [31] illustrated that QLV with a Mooney–Rivlin elastic strain energy function can fit data at a range of strain levels with relatively little error. Pioletti and Rakotomanana [32] quantified the ability of a time/strain separable model (such as QLV) to fit data over a range of tendon and ligament types and concluded that QLV was valid below certain specified strain levels that depend on the tissue type. Limbert and Middleton developed a finite strain viscoelastic constitutive law in which the elastic and viscous potentials were decoupled and used it to model human anterior cruciate [33] and posterior cruciate [34] ligaments. Johnson et al. [35] developed a finite strain phenomenological viscoelastic model for tendons and ligaments. Decreamer et al. [36] developed a microstructural viscoelastic model for the macroscopic behavior of soft tissue but the upscaling technique gave rise to QLV behavior. Other relevant papers include [3741], which are discussed in the context of this work in Sec. 4. For further references, we refer the reader to the reviews of tendon modeling by Reese and Weiss [42] and Thompson et al. [43].

An issue often raised with phenomenological models is that their parameters are not necessarily physical, whereas microstructural models often suffer from having too many parameters. The aim of this paper is to develop a microstructural model that relies on only a small number of parameters, all of which have a direct physical interpretation. Lanir [44] developed two models based on this principle: a high-density crosslink model and a low density crosslink (LDCL) model. In the high-density crosslink model, it was assumed that each viscoelastic collagen fibril is attached to a purely elastic elastin fiber, whereas in the LDCL model it was assumed that each viscoelastic collagen fibril moved independently through the extracollagenous matrix, which was assumed to be purely elastic. Sverdlik and Lanir [45] extended this model to three dimensions, using QLV for a single fibril and also incorporated the effects of preconditioning, and Raz and Lanir [46] further developed the model to eliminate the possibility of negative stresses in the tendon being predicted upon unloading. This approach predicts the results of certain experiments well, but does not predict strain-dependent relaxation and cannot be used to track the fibril length distribution in a tissue over time. In this paper, we use a model similar to the LDCL model and the model developed by Raz and Lanir, but here we also model the matrix as viscoelastic and use a more physically realistic method of monitoring the stress and strain in each collagen fibril, which we shall describe in detail in Sec. 2.3, by accounting for fibril creep. Due to the multiscale nature of this model, it is able to account for the stress relaxation that has been observed to take place at both the fibril and fiber levels [47,48]. The inclusion of the matrix phase gives rise to strain-dependent relaxation, and the inclusion of fibril creep allows us to predict the evolution of a fascicle's fibril length distribution over time for a given imposed fascicle strain.

The maximum level of strain induced in a tendon is relatively small (usually in the range 3–10%) and therefore, as a first approximation, a small strain linear model for its individual constituent phases is employed. By using a small strain model, we are implicitly neglecting $O(e2)$ terms, where e is the imposed strain, which would lead to a worst-case relative error of O(e) (i.e., 10% relative error at 10% applied strain), although as discussed later, since most fibrils within a tendon experience a strain lower than the macroscale strain, the actual relative error would be lower than this for any realistic initial fibril length distribution. During deformation, the fibrils interact with each other and the matrix as a result of their relative sliding [49]. It is assumed that this interaction induces an effective viscosity that results in the viscoelastic constitutive behavior of each component of the model. The structure of the paper is as follows. In Sec. 2, the model is introduced and the constitutive behaviors of the individual components are described, with particular attention being paid to their viscoelasticity. We then describe how the overall response of a fascicle is determined and consider two types of deformation: first, an incremental relaxation test and second, a cycle test. In Sec. 3, we describe experiments that were carried out to test the model and show that a good fit for both deformation types can be achieved with a single set of realistic parameter values for each fascicle. The fibril crimp distribution parameters of each fascicle are predicted by fitting them to the relaxation data, and then these parameters are used to predict the response to cyclic loading. We then demonstrate that the model is able to predict the time-evolution of the fibril length distribution for a given fascicle under cyclic loading. Finally, we show that the model predicts strain-dependent relaxation and gives a microstructrual explanation for the origin of this phenomenon. A discussion and conclusions are provided in Sec. 4.

## 2 The model

### 2.1 Configuration and Parameters.

We model the fascicle as a two-phase inhomogeneous medium consisting of a viscoelastic extracollagenous matrix phase and a large number of viscoelastic fibrils, each of which is initially crimped in the rest state of the fascicle. Our approach is based on the SSL model [10], the LDCL model described by Lanir [44], and the model of Raz and Lanir [46], in which the macroscopic response of a fascicle is due to the micromechanical behavior of the individual fibrils within it, which have varying lengths that are accommodated within the length of the fascicle via their crimp. It is assumed that the fibrils are coaligned with the fascicle and that they do not rotate. As the fascicle is stretched, the shortest fibrils are the first to become taut and start contributing to its stiffness. As the stretching continues, more fibrils straighten, and the stiffness of the fascicle increases, accounting for its exponential-shaped stress–strain curve. In our model, we additionally take into account the contribution to the macroscopic stress of the matrix phase. The fibrils are assumed to reside within the matrix and it is assumed that their interactions with each other and with the matrix cause the effective viscosity in each phase. We write the total volume fraction of (crimped and uncrimped) fibrils as $ϕ$ so that the volume fraction of the matrix phase is $1−ϕ$. In addition to the fibril volume fraction, we also need to know the relaxation behavior of the fibrils and the matrix, which shall be formally defined in terms of their relaxation functions later. The final input to the model is the distribution of fibril crimp lengths.

We model both the fibrils and the matrix as linear viscoelastic. We note that, as discussed above, linear viscoelasticity is inappropriate to model the averaged properties of tendon fascicles; however, since our model assumes a distribution of different fibril lengths within the fascicle, several nonlinear viscoelastic features arise as a result of the averaged behavior of these individual subunits.

In the following, we define the fibril and matrix stress and strain as $σf,ef$ and $σm,em$, respectively, which are scalars rather than tensors since this is a one-dimensional model. We denote the imposed macroscopic (fascicle) strain as e and the resulting macroscopic (fascicle) stress as σ. We determine the latter from what happens on the microscale. In other words, we formulate expressions relating the macroscopic stress to the fibril and matrix stresses at each point in time t. Note that the macroscopic strain e(t) is imposed for all time.

### 2.2 Matrix Response.

Since it is the simpler of the two phases, we begin by considering the response of the extracollagenous matrix to an applied fascicle load. We assume that the matrix strain is equal to the applied fascicle strain
$em(t)=e(t)$
(1)
and we model the matrix as a linear viscoelastic material, so that we have
$σm(t)=∫−∞tEm(t−s)demdsds$
(2)

$=Em(0)em(t)+∫0tEm′(t−s)em(s)ds$
(3)
where $Em(t)$ is the matrix relaxation function and $′$ represents differentiation with respect to the argument. The form (3), which is obtained from Eq. (2) upon using integration by parts, assumes the deformation begins at t =0 and is valid for all $t≥0$. We shall choose $Em(t)$ to take the form of a one-term Prony series
$Em(t)=E∞m+(E0m−E∞m)e−t/τrm$
(4)

where $E0m$ and $E∞m$ can be identified as the instantaneous and long-time Young's moduli of the matrix, respectively, and $τrm$ is its relaxation time.

### 2.3 Single Fibril Response.

Next, we consider the response of a single fibril to an applied fascicle strain. The length of this fibril when the fascicle is in its undeformed state shall be denoted $ℓ0$ and the associated undeformed fascicle length will be denoted L0 ($<ℓ0$), as depicted in Fig. 2. As the fascicle is stretched, the fibril length will change as a function of time, so that we can write the fascicle length at time t as $L(t)=L0+ΔL(t)$, and the corresponding fibril length as $ℓ(t)=ℓ0+Δℓ(t)$. The strains experienced by the fascicle and the fibril at time t can be defined as
$e(t)=L(t)−L0L0, ef(t)=ℓ(t)−ℓ0ℓ0$
(5)
respectively. Until taut, the fibril is passive and does not take up any stress; the only effect of the imposed displacement is to straighten out its crimp. At some time t = tA, however, the fibril will fully straighten, at which point we have $L(tA)=ℓ0$. The strain in the fascicle at this time is
$e(tA)=ec=ℓ0−L0L0$
(6)

We shall call this strain the critical strain, which we note is dependent on the length $ℓ0$ and so is unique to each individual fibril. We assume that the fibril has no resistance to compression along its longitudinal axis. If the fascicle is unloaded, the fibril will slacken as soon as $σf=0$ and therefore will never experience any negative stress. Due to viscoelastic memory effects, the critical fascicle strain at which straightening occurs changes on each subsequent reload.

From the moment that the fibril is first loaded until it slackens, we have $ℓ(t)=L(t)$ and so we can write the strain in the taut fibril as
$eft(t)=L(t)−ℓ0ℓ0=L0ℓ0(L(t)−L0L0−ℓ0−L0L0)=11+ec(e(t)−ec)$
(7)
where the superscript t is used to emphasize that the fibril is taut. As with the extracollagenous matrix, the fibril will be modeled as linear viscoelastic with relaxation function E(t), which we shall choose to take the form of a two-term Prony series, which has been shown to give a good fit to the relaxation behavior of individual collagen fibrils [3,5],
$E(t)=E∞+E1e−t/τr1+E2e−t/τr2$
(8)
in which $E∞$, E1, and E2 are relaxation moduli ($E∞$ is the long-time fibril Young's modulus and $E0=E∞+E1+E2$ is the instantaneous fibril Young's modulus) and $τr1$ and $τr2$ are relaxation times associated with the different relaxation time-scales, respectively. The stress in the taut fibril, therefore, is given by
$σft(t)=E(0)eft(t)+∫0tE′(t−s)ef(s)ds$
(9)

We note that the integral in this equation is over the entire strain history (assuming that the first fascicle deformation begins at t =0), not just over those times when the fibril is taut. This point is important for subsequent loadings.

If the fascicle is unloaded to the point that the fibril slackens, it will be important to monitor the strain within that fibril as it creeps as a function of time. To do this, we can invert the constitutive expression (9) to obtain
$ef(t)=J(0)σf(t)+∫0tJ′(t−s)σf(s)ds$
(10)
where J(t) is the fibril creep function, taking the form
$J(t)=J0+J1e−t/τc1+J2e−t/τc2$
(11)

where here, J0, J1, and J2 are creep moduli and $τc1$ and $τc2$ are creep times. The viscoelastic creep moduli and creep times can be expressed explicitly in terms of the relaxation moduli and times (see Appendix).

The strain in the slack fibril, therefore, is given by
$efs(t)=∫0tJ′(t−s)σf(s)ds$
(12)

where we have used the fact that $σf(t)=0$ while the fibril is slack and note that the integral is over the entire stress history. The superscript s emphasizes that the fibril is slack.

One must be aware of what is being imposed and what must be derived in this model. While taut, the fibril strain is imposed according to Eq. (7) and its stress must be derived from Eq. (9). This remains the case until the fibril stress reaches zero (as a result of reducing the fascicle strain), at which point the fibril slackens. After this point, the condition of zero fibril stress ($σf=0$) is imposed and its strain must be determined from Eq. (12). This will then remain the case until the fibril is reloaded to the point that it once again straightens out. At this point, the fibril strain is once again imposed and the fibril stress derived from Eq. (9). This is illustrated in Sec. 2.6 where we consider a cycle test consisting of two loading and unloading cycles.

At t =0, the fibril strain and stress are taken to be zero: $ef(0)=0, σf(0)=0$. For any subsequent time t, the fibril strain and stress are derived algorithmically via the flowchart in Fig. 3. This figure may be compared with Fig. 4, which shows the flowchart that would be used to derive the fibril stress and strain in the LDCL model derived by Lanir [44]. The functions in this flowchart are defined as follows:
$efL(t)=e(t)−ec, σfL=∫0tGL(t−s)H(efL(s))de(s)dsds$
(13)

where $GL(t)$ is the relaxation function associated with the LDCL model, and $H(·)$ is the unit step function. We note that Lanir correctly derived the fibril strain given by Eq. (7), but argued that this could be approximated by the expression in Eq. (13) since ec was assumed to be small. We also note that the LDCL model was derived in terms of forces, rather than stresses; however, since stress is simply force per unit area, the second expression in Eq. (13) is equivalent to the force–strain relationship of the LDCL model.

The flowchart in Fig. 4 can also be used to derive the fibril stress and strain in the model derived by Raz and Lanir [46], by instead using the following expressions:
$efL(t)=e(t)−ec1+2ec, σfL=∫0tGL(t−s)dσfedefLdefL(s)dsds$
(14)
where $σfe$ is the so-called “immediate” (or “elastic”) stress response, given by
$σfe={AefLefL≥00otherwise$
(15)

where A is the elastic coefficient of the fibrils. The factor of 2 in the equation for $efL$ is due to the fact that Raz and Lanir used finite Lagrangian strains rather than infinitesimal strains in their model.

By comparing these two flowcharts, the difference between the current model and the LDCL and Raz and Lanir models becomes clear. In the model described in this paper, the fibril strain is only imposed while the fibril is taut and the fibril creeps when it is unloaded, whereas in the LDCL and Raz and Lanir models, the fibril strain is imposed for all time and creep is unaccounted for. This can lead to negative fibril stresses in both the LDCL and the Raz and Lanir models when cyclic loading is considered.

### 2.4 Total Fascicle Response.

Once we know the stresses in the matrix and in a fibril that tautens at critical fascicle strain ec, it is simple to obtain the homogenized response of a fascicle containing fibrils with a distribution of crimp by volume averaging. The total fascicle stress is equal to the average of the stresses in the fibrils and matrix in proportion to their volume fractions; i.e., $(1−ϕ)$ times the matrix stress plus $ϕ$ times the average fibril stress, where $ϕ$ is the fibril volume fraction. Since each fibril has its own associated critical strain, the average fibril stress is determined by integrating over the individual fibril stresses multiplied by a critical strain distribution function. Considering a unit volume of fascicle we have
$σ(t)=(1−ϕ)σm(t)+ϕ∫0∞p(ec)σf(e(t),ec)dec, t∈[0,∞)$
(16)

where $p(ec)$ is the initial distribution of critical strains and the notation $σf(e(t),ec)$ is used to emphasize the fact that the stress in each individual fibril is dependent on both the imposed fascicle strain, e(t), and the critical fascicle strain, ec, at which that fibril first becomes taut.

We note that, in the model proposed here, any continuous function may be chosen for the imposed fascicle strain, provided it satisfies $e(0)=0$. Equation (16) provides a complete description of the model for any deformation, but the main complexity in its implementation arises in calculating the fibril stress σf according to the flowchart in Fig. 3. The specific form of the deformation being considered determines how challenging it is to calculate σf; therefore, to illustrate the implementation of the model, we consider two tests: a relaxation test and a cycle test.

### 2.5 Relaxation Test.

The first example that illustrates the model is a relaxation test in which the fascicle strain is rapidly ratcheted up in several discrete steps with a long period of constant strain in between each step, as depicted in Fig. 5(a). In Fig. 5(b), the strain in a fibril with associated critical strain $ec>0$ is displayed. The critical time tA at which the fibril first completely straightens is labeled. In Fig. 5(c), the corresponding fibril stress is shown. Since the fascicle is not unloaded in this problem, the fibril strain is 0 for $t≤tA$ and is defined by Eq. (7) for $t≥tA$ and the fibril stress is determined by Eq. (9) for all $t≥0$. The matrix strain and stress are defined by Eqs. (1) and (3), respectively, for all $t≥0$. The resulting fascicle stress (16) will be plotted and compared with experimental data in Sec. 3.2.

### 2.6 Cycle Test.

The second example is a cycle test, as depicted schematically in Fig. 6. Figure 6(a) illustrates the imposed fascicle strain, and Figs. 6(b) and 6(c) show the corresponding strain and stress, respectively, in a single fibril of length $ℓ0>L0$. We note that this test is more challenging to model than the relaxation test due to the unloading that occurs between tB and tD and $tB′$ and $tD′$. The matrix strain and stress are once again defined by Eqs. (1) and (3), respectively, for all $t≥0$; however, to calculate the fibril strain and stress, we must consider each stage of the cycle test separately as we now discuss.

The response from t =0 through t = tA to the maximum strain $e(tB)$ is due to the initial loading, a single fibril only starts to take up stress at time tA when it has fully straightened. At this point, $e(tA)=ec$, but $ef(tA)=0$. We see in Fig. 6(b) that on this initial load curve
$ef(t)=e(t)−ec1+ecH(e(t)−ec)$
(17)
where the unit step function ensures that the strain only begins to increase once the critical strain is reached and σf begins to increase after t = tA. Along this first load curve, the fibril stress is therefore calculated using Eq. (9) in the form
$σf(ef) ={0t∈[0,tA]E(0)ef(t)+∫0tE′(t−s)ef(s)dst∈[tA,tB]$
(18)
although it is more useful to use the expression (17) and write
$σf(e(t),ec) ={0t∈[0,tA]E(0)e(t)−ec1+ec+∫tAtE′(t−s)e(s)−ec1+ecdst∈[tA,tB]$
(19)

Since ec and e(t) are known, the fibril stress is fully specified in any fibril up to t = tB.

Once the maximum strain is achieved at t = tB, the macroscopic strain is then reduced at a specified strain rate back to zero. Initially, the fibril strain will continue to be defined by Eq. (7); however, at some time t = tC, critically, the fibril stress becomes zero. We will discuss this shortly, but first we state, therefore, that the range of times for which (19) applies is extended so that
$σf(e(t),ec) ={0t∈[0,tA]E(0)e(t)−ec1+ec+∫tAtE′(t−s)e(s)−ec1+ecdst∈[tA,tC]$
(20)
Along the first part of the unload curve (from t = tB to t = tC), the fibril stress will not be the same as the stress for the corresponding stretch ef on the load curve due to stress relaxation and the memory effect. Its qualitative behavior in this first unloading curve is illustrated in Fig. 6(c). As a result of this effect, for each fibril, at some time t = tC the fibril stress will become zero. As previously discussed, our assumption is that the fibrils are unable to take up compressive stress and therefore, at this point the fibril will slacken, maintaining zero stress within the fibril, and the crimp will begin to reappear, so that as time progresses we are required to determine the fibril strain from the zero stress condition at time t. In other words, after tC (up until some, as yet unspecified, time, $tA′$ say) we must use (10) to determine the time evolution of the strain as a function of the stress, i.e.,
$ef(t)=∫tAtCJ′(t−s)σf(s)ds, t∈[tC,tA′]$
(21)
while imposing zero stress
$σf(e(t),ec)=0, t∈[tC,tA′]$
(22)

The upper limit of integration in Eq. (21) is tC since the fibril stress is zero after this time, but the fibril strain evolves in time, of course, due to the time dependence of the creep function in the integrand. This expression for the fibril strain is valid until some time $t=tA′$ as depicted on Fig. 6(a), as we will discuss shortly.

While the fibril stress is zero on the unload curve, the macroscopic strain e(t) is itself being reduced back to zero, i.e., to t = tD on the curve in 6(a). If we maintained e(t)=0 for all $t≥tD$, the fibril strain $ef(t)$ in Eq. (21) would also reduce back to zero as $t→∞$; however, we choose instead to immediately reload as indicated in the figure, so that the macroscopic strain begins to increase once again. At this time, the fibril stress is still zero, the fibril is crimped but it still has within it some nonzero strain ef. We now discuss the re-load curve and how we determine the time $tA′$ at which the stress attains a positive value once again.

After reducing the macroscopic strain to zero at t = tD, it is once again increased for the second load cycle. In this illustrative example, we have used the same strain rate for the second load cycle, but the theory would be valid for any other strain rate. At t = tD, the fibril strain is still reducing to accommodate the zero fibril stress condition; however, at some later time $t=tA′$ when
$e(tA′)=(1+ec)ef(tA′)+ec$
(23)

the fibril will start to take up stress once again since it will have straightened (this expression is obtained by substituting $t=tA′$ into Eq. (7) and rearranging). We note that this second critical strain is larger than the first. This is because the fibril had an initial “rest length” $ℓ0$, but, as a result of the first load cycle, it is still strained in its zero stress state. The “rest length,” therefore, has increased in its slackened state (we can call this new rest length $ℓ1(t)$, which is a monotonically decreasing function of time that tends to 0 as $t→∞$). During the reload curve, the fibril will straighten when the fascicle length L(t) is equal to the new rest length. The second critical strain, therefore, is larger than the first for all fibrils that underwent straightening in the first load cycle. For those fibrils that did not undergo straightening, the second critical strain is equal to the first.

$σf(e(t),ec)=E(0)e(t)−ec1+ec+∫tAtE′(t−s)ef(s)ds, t∈[tA′,tC′]$
(24)
which extends to $tC′$ for the same reason as the first loading. As such, in summary, over the full two-cycles the fibril stress is given by
$σf(e(t),ec)={0t∈[0,tA]E(0)e(t)−ec1+ec+∫tAtE′(t−s)ef(s)dst∈[tA,tC]0t∈[tC,tA′]E(0)e(t)−ec1+ec+∫tAtE′(t−s)ef(s)dst∈[tA′,tC′]0t∈[tC′,∞)$
(25)
and the fibril strain is given by
$ef(e(t),ec) ={0t∈[0,tA]e(t)−ec1+ect∈[tA,tC]∫tAtCJ′(t−s)σf(s)dst∈[tC,tA′]e(t)−ec1+ect∈[tA′,tC′]∫tAtCJ′(t−s)σf(s)ds+∫tA′tC′J′(t−s)σf(s)dst∈[tC′,∞)$
(26)

We can see from the above expressions that the stress and strain in a given fibril are dependent only on the macroscopic imposed fascicle strain e(t) and on the critical strain ec at which it first becomes taut. Each fibril will have a different secondary critical strain which will depend on the deformation history, but this can be determined directly from the initial critical strain and the imposed fascicle strain.

## 3 Experiments and Model Predictions

### 3.1 Materials and Methods

#### 3.1.1 Rat Tail Tendon Preparation.

Collagen fascicles were extracted from the tail of a 3-month-old female Sprague-Dawley rat, killed for other, nonrelated reasons, following previously published protocols [50]. The skin was resected to expose the four tail tendons, and fascicles extracted from the proximal end of a single superior quadrant. The fascicles were teased directly from their endotendinous sheaths, taking care not to load them during the dissection process, and cut transversely to lengths of approximately 60 mm. Fascicles were dissected immediately before use, to prevent drying prior to analysis.

#### 3.1.2 Mechanical Assays.

Fascicles were tested in previously described custom-made loading chambers, allowing full hydration of samples throughout testing [51]. The diameter of each fascicle was first determined, taking continuous readings along the 10 mm test length of the fascicle using a laser micrometer (LSM-501, Mitutoyo, Kawasaki, Japan), from which the smallest diameter was used to calculate the cross-sectional area, assuming a circular shape. The measured diameters were 178 μm for fascicle 1, 205 μm for fascicle 2, and 186 μm for fascicle 3. The fascicle was then secured between the stainless steel grips of the loading chamber at a test length of 10 mm, the chamber filled with Dulbecco's modified Eagle's medium and sealed, and the chamber secured within a material test machine loading frame (InstronElectroPuls1000, Instron, High Wycombe, UK) equipped with a 250 N load cell.

A preload of 0.02 N was applied to each fascicle to provide a consistent starting point for the test, and the new distance between the grips recorded as the effective gage length, with sample strains calculated accordingly. Samples were preconditioned with 20 triangular waveform cycles at 1 Hz between 0 and 6% strain, then held for 1 h at 0% strain to ensure all fibrils had fully relaxed prior to initiating the test, which consisted of ten loading cycles followed by an incremental stress relaxation test. The ten loading cycles were applied at 1 Hz between 0 and 6% strain (triangular waveform), after which the sample was held at 0% strain for an hour to allow recovery. The incremental stress relaxation test consisted of three strain increments, taking the sample to 2%, 4%, and 6% strain at a rate of 12% per second, holding the sample at each increment for 5 min, before progressing immediately to the next strain increment. Force and displacement data were collected at 100 Hz during testing, and stress-time plots produced for each sample.

### 3.2 Comparison Between Theory and Experiments

#### 3.2.1 Parameter Selection.

To calculate the theoretical stress predicted by the model, we must select values for the fibril and matrix constitutive parameters, and also to select a function for the initial critical strain distribution $p(ec)$. There are two studies that fit Prony series relaxation functions to relaxation tests on individual collagen fibrils—one tested sea cucumber collagen [3] and the other tested bovine Achilles collagen [5]. Since bovine collagen is mammalian (and therefore likely to have mechanical behavior that is more similar to rat tail tendon), the parameter values reported by Yang et al. [5] for native collagen fibrils were used in our model. In that paper, the relaxation moduli were reported in nondimensional form as follows:
$A1=E1E∞+E1+E2=0.2, A2=E2E∞+E1+E2=0.17$
(27)
Therefore, to complete our relaxation function, we still need to independently determine one of the three relaxation moduli. We therefore used the literature on rat tail tendon fibril Young's moduli to determine the value of $E∞$. The range of collagen Young's moduli for rat tail tendon fibrils reported in the literature is wide, ranging from tens of MPa [52] to 11.5 GPa [53]; however, the value chosen for this study was 3 GPa, which is the lower of the values measured by Grant et al. [54] in a recent study and gives a value of $E0=E∞+E1+E2$ that is close the upper value of 5 GPa report by Grant et al. Therefore, a summary of the fibril parameters used here is as follows (quantities reported to two significant figures):
$E∞=3.0 GPa, E1=950 MPa, E2=810 MPaτr1=1.9 s, τr2=52 s$
(28)
Note that these parameters were selected a priori and were not used as fitting parameters. To the authors' knowledge, there is no published study that reports the Young's modulus of tendon extracollagenous matrix explicitly; however, Henninger et al. [55] reported that the elastin (one of the major components of tendon matrix) contribution to the peak stress of porcine medial collateral ligaments under 10% axial strain is approximately 2 MPa (implying an effective Young's modulus of approximately 20 MPa). The matrix shear modulus has been estimated to be even lower, at approximately 3–5 kPa [56]. Since these values are so low compared to the fibril moduli, we have chosen to neglect the matrix contribution in our initial modeling of the relaxation and cycle tests. This decision reduces computational time and the number of parameters in the model without significantly affecting the final reported fascicle stresses. We later re-introduce the matrix phase in order to explain strain-dependent relaxation. The collagen volume fraction was selected to be
$ϕ=0.80$
(29)

This was the collagen area fraction (to two significant figures) in 4–6-month-old Wistar rat tail tendons determined by Screen et al. [57]; we assume the collagen volume fraction in rat tail tendon does not differ significantly from this value.

Finally, the initial critical strain distribution was modeled by a truncated normal distribution, that it is defined over the range $(0,∞)$
$p(ec)=e−(ec−μ)22s22πs(1−12erfc(μ2s))$
(30)

where μ and s are the truncated equivalents of the mean and standard deviation of the normal probability distribution function, respectively. This distribution was chosen due to the fact that its parameters can be understood intuitively. The values of μ and s were predicted by fitting them to the relaxation test data, as discussed below.

#### 3.2.2 Relaxation Tests.

To model the incremental relaxation tests, the following fascicle strain was input into the model:
$e(t)={0.12t0≤t<0.170.020.17≤t<300.170.02+0.12(t−300.17)300.17≤t<300.330.04300.33≤t<600.330.04+0.12(t−600.33)600.33≤t<600.50.06600.5≤t<900.5$
(31)

where t is measured in seconds.

To remove noise and obtain data suitable for fitting, a 1000 time-step (equivalent to 10 s) moving average was applied to each dataset. The values of μ and s were determined by minimizing the l2 norm of the residuals between the theoretical and experimental values at three points along the averaged relaxation curves, as illustrated in Fig. 7.

#### 3.2.3 Cycle Tests.

For the cycle tests, the first two cycles of the ten that were carried out were modeled. The exact form for the fascicle strain used in the model was
$e(t)={0.12t0≤t<0.50.06−0.12t0.5≤t<10.12(t−1)1≤t<1.50.06−0.12(t−1)1.5≤t<2$
(32)

where again, t is measured in seconds.

For the cycle tests, the values of μ and s obtained from the relaxation tests were used. Therefore, these theoretical plots are predictions as no fitting took place to match the cycle test results directly.

#### 3.2.4 Results.

The values of μ and s used to fit each relaxation test are shown in Table 1. The results of the fitting to the relaxation tests are shown in Figs. 8(a)8(c) (these plots show the raw data, as opposed to the averaged data shown in Fig. 7), and the resulting predictions for the cycle tests are shown in Figs. 8(d)8(f). We note the remarkable agreement with the cyclic data given that these datasets did not contribute to the fitting procedure. In cycle tests 2 and 3, a small amount of friction was evident at the bottom of the actuator travel, where the top grip came into contact with the chamber; therefore, in order to minimize its effects on the data, a graded moving average filter (over nine time points close to the peaks, increasing gradually to 15 time points in the troughs) was used to smooth the stress data at the bases of the troughs for all three datasets. The predicted critical strain distribution functions are shown in Fig. 9.

To assess the precision of the model's predictions, three contour plots are provided in Fig. 10 which show the sum of the squared residuals at the three points illustrated in Fig. 7. As can be seen, the range of values of μ and s that fit the first relaxation test is relatively wide compared to those for the second and third relaxation tests.

### 3.3 Time-Evolution of Fibril Length Distribution.

Since the model presented in this paper accounts for fibril creep, for a given imposed fascicle strain e(t) and a given initial fibril length 0, we can calculate the associated fibril strain $ef(t)$, and therefore, the current fibril length l(t) at any point in time. An animation of the time-evolution of the lengths of 1000 fibrils randomly sampled from a fascicle with $μ=0.07$ and s =0.03 that is subjected to the deformation described in Eq. (32) can be found in the data given in the footnote link along with an animation of the case when fibril creep is neglected.2 The length distribution is plotted in nondimensional form as fibril length divided by initial fascicle length (i.e., $l(t)/L$).

### 3.4 Strain-Dependent Relaxation.

A key benefit of the model described above is that it predicts that the relaxation rate of a tendon depends upon the magnitude of the strain—a prediction that is backed up by experimental data [22,23]. To demonstrate this prediction, however, it is necessary to re-introduce the matrix phase. In the absence of any experimental mechanical data on the viscoelastic behavior of this phase, we shall assume that its instantaneous Young's modulus is equal to 20 MPa, as derived previously from the work of Henninger et al. [55], and that the long-term modulus is half of that value at 10 MPa. We shall also assume that the relaxation time of the matrix is much shorter than that of the fibrils—we assume it to be one second. A summary of these parameter values is as follows:
$E0m=20 MPa, E∞m=10 MPa, τrm=1 s$
(33)
Using these parameters, along with those for the fibrils in Eq. (28), we plot the relaxation behavior of a fascicle with $μ=0.05$, s =0.004 that has been stretched to three different initial strains. The assumed deformation is
$e(t)={emax(t/0.01)0≤t<0.01emaxt≥0.01$
(34)

where $emax$ is the maximum strain. Note the rapid initial strain rate. We consider three different values of $emax$: $2%, 4%$, and 6%, and in Fig. 11 we plot the corresponding nondimensionalized stresses $σ(t)/σ(0.01)$ as dashed lines. The solid lines for comparison are taken from the experimental data plotted in Fig. 7, by taking the mean of the relaxation responses of the three fascicles at each initial strain level and nondimensionalizing the data on the peak stress. A 100 time-step moving average was applied over the first 500 data points, and a 1000 time-step moving average was applied to the rest of the data in order to remove noise. As can be seen, the model qualitatively predicts the behavior that is observed experimentally. We emphasize that the LDCL and Raz and Lanir models do not predict this behavior—the nondimensionalized relaxation curves would be identical for all three initial strains for both of these models.

## 4 Discussion

We have developed a model of tendon viscoelasticity based on fascicle microstructure. The fibrils and extracollagenous matrix were assumed to be linear viscoelastic, and it was shown that the complex nonlinear behavior exhibited by tendons can be explained as a geometrical effect. In other words, the findings of the model are consistent with the hypothesis that tendon fibrils and matrix are fundamental units that always have the same mechanical properties (within a given tendon, at least) and that differences between the behaviors of different fascicles are caused by differences in the distributions of the lengths of their fibrils.

The model agrees well with experiments and is able to fit multiple datasets with a single set of fibril constitutive parameters. Additionally, we note that the results plotted in this paper do not necessarily represent the best fit that could be achieved with this model. Given the wide range of reported values for the fibril Young's modulus in the literature, it is impossible to be certain that the values used in this paper are correct. Instead, the parameters used represent an educated choice given the information available. We emphasize they were decided upon a priori and were not used as fitting parameters. It is likely that if it were possible to determine the true values of these parameters with certainty, then the model would fit the mechanical test data even better. To achieve this aim, it would be useful to carry out relaxation tests on individual rat tail tendon collagen fibrils similar to those performed on sea cucumber fibrils by Shen et al. [3] and on porcine Achilles fibrils by Yang et al. [5]. Nevertheless, the fact that the complete set of parameters (including the distribution parameters that were determined by fitting to the relaxation data) can be used to predict the outcome of the cycle tests so accurately (without requiring any fitting to the cycle test data itself) demonstrates the predictive power of this model. We note that both the peak stresses and the gradients of the curves are predicted accurately. Thus, it appears that the model captures the essential physics that governs the viscoelastic behavior of tendon fascicles.

The fibril critical strain distributions used to fit the data in this paper should not be seen as predictions of their true shapes. It would only be possible to make such predictions if we had a much higher degree of certainty in the fibril parameters. An increase in the fibril relaxation moduli can be compensated for by moving the mean of the truncated normal distribution further from zero and/or by modifying the standard deviation while still achieving a good fit to the data. Conversely, a decrease in the fibril relaxation moduli can be compensated for by bringing the mean closer to zero. The shape of the predicted critical strain distribution function for fascicle 1 appears to be qualitatively different to those of fascicles 2 and 3; however, it is interesting to note that Fig. 10(a) shows that the range of values of μ and s that fit this dataset well is relatively wide. By following the yellow contour, it can be seen that it is still possible to fit the data well with a value of μ much closer to those predicted for fascicles 2 and 3. If a value closer to $μ=0.04$ was imposed, for example, the distribution function for fascicle 1 would look much more similar to those of fascicles 2 and 3. Another possibility is that the true fibril critical strain distributions are bi- or multimodal. We did not consider these cases as we wanted to keep the number of fitting parameters in our model to a minimum; however, if these possibilities were considered, it would be possible to improve the fit to the data significantly.

While our model and datasets do not allow us to make definitive predictions about the values of the fibril constitutive parameters, they do allow us to put a lower bound on one of their values. For fascicle 2, the jump in stress between the 4% increment and 6% strain increment once fully relaxed is approximately 36 MPa. Even if all of the fibrils were taut by this point and the collagen volume fraction was 1, such a jump would require a long-time fibril Young's modulus of $E∞=36/0.02 MPa=1800$ MPa. Therefore, we can state with confidence that the true value of this parameter is strictly greater than 1800 MPa, provided that the experimental measurements and the model are accurate. One potential source of error in stress measurements of a tendon fascicle is the fact that its cross-sectional area varies along its length. In our experiments, the smallest measured diameter was used to derive the cross-sectional area assuming a circular cross section. Assuming a circular cross section gives an overestimate of the true cross-sectional area of approximately 4% [58], which is why the smallest measured diameter was used.

The inclusion of fibril creep in our model allowed us to predict the time-evolution of the fibril length distribution in a fascicle. The ability to determine individual fibril strain histories may be important for predicting tendon rupture under a given time-dependent deformation; therefore, fibril creep is a crucial feature of the model that offers an improvement upon previous viscoelastic collagen recruitment models. As can be seen in the data given in the footnote link,2 neglecting fibril creep results in a significant number of the fibrils being compressed to be shorter than their initial length during the unloading phase of a cyclic test, which would give rise to negative fibril stresses. This would not be expected in crimped fibrils, and the resulting length distribution functions have an unusual shape when fibril creep is neglected. Bevan et al. [59] recently measured the recruitment probability density function of tendon fibrils under quasi-static deformation conditions using confocal microscopy. They used a triangular distribution to fit their data and predicted a modal critical strain of around 1% and a maximal critical strain of around 3%. These predicted values are slightly smaller than those predicted on average in this paper, but are of the same order of magnitude. If it were possible to adapt this method to arbitrary strain rates, that could be a way to validate the fibril length evolution predictions of the proposed model.

As mentioned previously, one of the benefits of the proposed model is that it can account for the strain-dependent relaxation that has been observed experimentally in the literature [22,23]. Although other constitutive models have predicted strain-dependent relaxation (e.g. [60]), our approach has allowed us to provide a microstructural explanation for the origin of this phenomenon. It is interesting to note that the model predicts that this behavior is due to the difference in the relaxation times between the matrix and the fibrils. If we omit the matrix phase, as in Sec. 3.2, then this behavior is not exhibited, as indeed it is not in the LDCL and Raz and Lanir models. When the matrix is included, however, the relaxation behavior of the fascicle is similar to that of the matrix for small strains (hence the short relaxation time) and then becomes more and more like that of the fibrils as the initial strain is increased (hence the slower relaxation for higher strains). This may give some insight into the origins of strain-dependent relaxation behavior in other composite materials. Given the material parameters used here, the extent of relaxation is greater for small strains due to the difference between the initial and long-time Young's moduli of the matrix being bigger than that of the fibrils. We note, however, that, in the absence of experimental relaxation data on isolated extracollagenous matrix, we chose the matrix long-time Young's modulus and the relaxation time somewhat arbitrarily in order to illustrate the phenomenon of strain-dependent relaxation. Given these parameter values, the matrix relaxes more quickly than the fibrils, and therefore at small strains, the fascicles relax more quickly than at larger strains. If the matrix parameters were such that it relaxed more slowly than the fibrils, then this trend would reverse. Similarly, by changing the ratio of $E0m$ to $E∞m$, the extent of relaxation at small strains could be changed to be greater or lesser than that at large strains. If the nondimensional relaxation function of the matrix were identical to that of the fibrils, then the relaxation would no longer be strain-dependent.

In addition to the fibril and matrix constitutive parameters, in the long term it will be necessary to independently measure the geometrical parameters $ϕ$, μ, and s. These could potentially be measured by confocal microscopy [59,61] or X-ray computed tomography [6264]. If it were possible to measure all of the constitutive and geometrical parameters independently for a given fascicle that had also been mechanically tested, then the model could be fully validated.

There have been several other approaches to modeling the viscoelastic behavior of ligaments and tendons. De Vita and Slaughter [40] developed a structural model of the strain rate-dependent behavior of anterior cruciate ligaments that is valid for large strain and Sopakayang and De Vita [41] developed a model that is capable of accounting for creep, relaxation and strain stiffening in soft tissues. The formulation of the latter model was similar to that presented here and gave a good fit to several datasets. One key difference in that model is that the fibrils were modeled as elastic and the matrix was modeled via a Maxwell model. Given the evidence of the viscoelasticity of individual fibrils [3,5], however, we considered it important to include this feature within our model. In an early paper, Thornton et al. [38] showed that linear viscoelasticity cannot explain the creep behavior of rabbit medial collateral ligaments and the same group later went on to develop a structural model of ligament creep based on fibril recruitment [39]. This model was based on similar principles to that derived here; however, since the focus was entirely on creep, and they did not consider relaxation or cycle tests, they avoided the need to consider the transition of a fibril between relaxation and creep, which is one of the key contributions of the model presented in this paper (as illustrated schematically in Fig. 3). Together, the results of Thornton et al. along with those presented here indicate that fibril recruitment may play a role in both creep and relaxation.

A possible extension to our model would be to include more phases than the two considered here. In reality, the extracollagenous matrix is not a single phase, but is several materials (such as elastin, proteoglycans and tenocytes, for example,). It may be of interest in the future to model the interactions between these materials explicitly rather than bundling them together into a matrix phase. Another possible extension would be to model the interaction between the fibrils and the matrix. Ciarletta and Ben Amar [37] developed a dissipative theory for the bridges that arise between collagenous and noncollagenous proteins in the extracellular matrix, using an exponential strain energy function to account for the shape of the tissue's stress–strain curve phenomenologically. An interesting challenge would be to combine this theory with the microstructural viscoelastic framework presented in this paper to produce a model that accounts for both fibril–matrix interaction and fibril length distributions.

We conclude by noting that, although this model was developed for ligaments and tendons, it could be extended to higher dimensions in order to model other soft tissues. This could be achieved by introducing a probability distribution function to account for the alignment of the fibrils in addition to their critical strain distribution function. It could also be extended to model polymers by assuming that the fundamental unit in the model is a polymer chain oriented in three-dimensional space. For this to work, however, it would be necessary to track the strain in each individual polymer chain, which is likely to be computationally expensive.

## Acknowledgment

The authors are grateful to Dr. Riccardo De Pascalis who took part in some early discussions associated with this research and to Dr. Jean-Marc Allain who carried out experiments for an early draft of this paper. TS and WJP would like to thank the EPSRC for funding this work through fellowship grants EP/L017997/1 and EP/L018039/1, respectively. BL acknowledges the School of Mathematics at the University of Manchester (via the MAPLE Platform Grant (EP/I01912X/1)) for funding her research visit to the UK in April 2014. IDA is grateful to the Royal Society for funding his Wolfson Research Merit award (2013-2018), and to EPSRC/UKRI for grant EP/R014604/1.

## Funding Data

• Engineering and Physical Sciences Research Council (Grant Nos. EP/I01912X/1, EP/L017997/1, EP/L018039/1, and EP/R014604/1; Funder ID: 10.13039/501100000266).

• Royal Society (Wolfson Research Merit Award; Funder ID: 10.13039/501100000288).

## Footnotes

2

The raw experimental data is available to download and has the https://doi.org/10.15127/1.308061. Additional Supplementary Material has the https://doi.org/10.17632/yzzczmsttx.3. This includes the Mathematica codes used to fit and predict the experimental data and the animations referred to in Sec. 3.3.

### Creep Moduli and Relaxation Times

The creep moduli and relaxation times can be expressed in terms of the relaxation moduli and relaxation times by taking the Laplace transform of the constitutive equation, rearranging, then taking the inverse Laplace transform and comparing coefficients. The resulting relationships are as follows:
$J0=1E∞$
(A1)

$J1=E∞(E1−E2)(τr1−τr2)+(E1+E2)(E1τr1+E2τr2−α)2E∞E0(E1τr1+E2τr2+E∞(τr1+τr2))2−4E∞E0τr1τr2$
(A2)

$J2=E∞(E1−E2)(τr1−τr2)+(E1+E2)(E1τr1+E2τr2+α)2E∞E0(E1τr1+E2τr2+E∞(τr1+τr2))2−4E∞E0τr1τr2$
(A3)

$τc1=2E0τr1τr2(E∞+E1)τr1+(E∞+E2)τr2+α$
(A4)

$τc2=2E0τr1τr2(E∞+E1)τr1+(E∞+E2)τr2−α$
(A5)
where
$α=[E∞2(τr1−τr2)2+2E0(τr1−τr2)(E1τr1−E2τr2)+(E1τr1+E2τr2)2]12$
(A6)

## References

References
1.
Lin
,
T. W.
,
Cardenas
,
L.
, and
Soslowsky
,
L. J.
,
2004
, “
Biomechanics of Tendon Injury and Repair
,”
J. Biotech.
,
37
(
6
), pp.
865
877
.10.1016/j.jbiomech.2003.11.005
2.
Kastelic
,
J.
,
Galeski
,
A.
, and
Baer
,
E.
,
1978
, “
The Multicomposite Structure of Tendon
,”
Connect. Tissue Res.
,
6
(
1
), pp.
11
23
.10.3109/03008207809152283
3.
Shen
,
Z. L.
,
Kahn
,
H.
,
Ballarini
,
R.
, and
Eppell
,
S. J.
,
2011
, “
Viscoelastic Properties of Isolated Collagen Fibrils
,”
Biophys. J.
,
100
(
12
), pp.
3008
3015
.10.1016/j.bpj.2011.04.052
4.
Sopakayang
,
R.
,
De Vita
,
R.
,
Kwansa
,
A.
, and
Freeman
,
J. W.
,
2012
, “
Elastic and Viscoelastic Properties of a Type I Collagen Fiber
,”
J. Theor. Biol.
,
293
, pp.
197
205
.10.1016/j.jtbi.2011.10.018
5.
Yang
,
L.
,
van der Werf
,
K. O.
,
Dijkstra
,
P. J.
,
Feijen
,
J.
, and
Bennink
,
M. L.
,
2012
, “
Micromechanical Analysis of Native and Cross-Linked Cololagen Type I Fibrils Supports the Existence of Microfibrils
,”
J. Mech. Behav. Biomed. Mater.
,
6
, pp.
148
158
.10.1016/j.jmbbm.2011.11.008
6.
Hukins
,
D. W. L.
,
Leahy
,
J. C.
, and
Mathias
,
K. J.
,
1999
, “
Biomaterials: Defining the Mechanical Properties of Natural Tissues and Selection of Replacement Materials
,”
J. Mater. Chem.
,
9
(
3
), pp.
629
636
.10.1039/a807411i
7.
Rigby
,
B.
,
Hirai
,
N.
,
Spikes
,
J.
, and
Eyrinf
,
H.
,
1959
, “
The Mechanical Properties of Rat Tail Tendon
,”
J. Gen. Physiol.
,
43
(
2
), pp.
265
283
.10.1085/jgp.43.2.265
8.
Fung
,
Y. C. B.
,
1967
, “
Elasticity of Soft Tissues in Simple Elongation
,”
Am. J. Physiol.
,
213
(
6
), pp.
1532
1544
.10.1152/ajplegacy.1967.213.6.1532
9.
Holzapfel
,
G. A.
,
Gasser
,
T. C.
, and
Ogden
,
R. W.
,
2000
, “
A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models
,”
J. Elast.
,
61
(
1/3
), pp.
1
48
.10.1023/A:1010835316564
10.
Kastelic
,
J.
,
Palley
,
I.
, and
Baer
,
E.
,
1980
, “
A Structural Mechanical Model for Tendon Crimping
,”
J. Biotech.
,
13
(
10
), pp.
887
893
.10.1016/0021-9290(80)90177-3
11.
Lanir
,
Y.
,
1983
, “
Structure-Strength Relations in Mammalian Tendon
,”
J. Biotech.
,
16
(
1
), pp.
1
12
.10.1016/0021-9290(83)90041-6
12.
Ault
,
H. K.
, and
Hoffman
,
A. H.
,
1992
, “
A Composite Micromechanical Model for Connective Tissues: Part I–Theory
,”
ASME J. Biotech. Eng.
,
114
(
1
), pp.
137
141
.10.1115/1.2895437
13.
Ault
,
H. K.
, and
Hoffman
,
A. H.
,
1992
, “
A Composite Micromechanical Model for Connective Tissues—Part II: Application to Rat Tail Tendon and Joint Capsule
,”
ASME J. Biotech. Eng.
,
114
(
1
), pp.
142
146
.10.1115/1.2895438
14.
Shearer
,
T.
,
2015
, “
A New Strain Energy Function for the Hyperelastic Modelling of Ligaments and Tendons Based on Fascicle Microstructure
,”
J. Biotech.
,
48
(
2
), pp.
290
297
.10.1016/j.jbiomech.2014.11.031
15.
Beskos
,
D.
, and
Jenkins
,
J.
,
1975
, “
A Mechanical Model for Mammalian Tendon
,”
ASME J. Appl. Mech.
,
42
(
4
), pp.
755
758
.10.1115/1.3423699
16.
Freed
,
A.
, and
Doehring
,
T.
,
2005
, “
Elastic Model for Crimped Collagen Fibrils
,”
ASME J. Biotech. Eng.
,
127
(
4
), pp.
587
593
.10.1115/1.1934145
17.
Grytz
,
R.
, and
Meschke
,
G.
,
2009
, “
Constitutive Modeling of Crimped Collagen Fibrils in Soft Tissues
,”
J. Mech. Behav. Biotech.
,
2
(
5
), pp.
522
533
.10.1016/j.jmbbm.2008.12.009
18.
Shearer
,
T.
,
2015
, “
A New Strain Energy Function for Modelling Ligaments and Tendons Whose Fascicles Have a Helical Arrangement of Fibrils
,”
J. Biotech.
,
48
(
12
), pp.
3017
3025
.10.1016/j.jbiomech.2015.07.032
19.
Hurschler
,
C.
,
Loitz-Ramage
,
B.
, and
Vanderby
,
C.
,
1997
, “
A Structurally Based Stress-Stretch Relationship for Tendon and Ligament
,”
ASME J. Biotech. Eng.
,
119
(
4
), pp.
392
399
.10.1115/1.2798284
20.
Viidik
,
A.
,
1968
, “
A Rheological Model for Uncalcified Parallel-Fibered Collagenous Tissue
,”
J. Biotech.
,
1
(
1
), pp.
3
11
.10.1016/0021-9290(68)90032-8
21.
Viidik
,
A.
,
1973
, “
Functional Properties of Collagenous Tissues
,”
Int. Rev. Connect Tissue Res.
,
6
, pp.
127
215
.10.1016/B978-0-12-363706-2.50010-6
22.
Provezano
,
P.
,
Lakes
,
R.
,
Keenan
,
T.
, and
Vanderby
,
R.
, Jr.
,
2001
, “
Nonlinear Ligament Viscoelasticity
,”
Ann. Biomed. Eng.
,
29
, pp.
908
914
.10.1114/1.1408926
23.
Duenwald
,
S. E.
,
Vanderby
,
R.
, Jr.
, and
Lakes
,
R. S.
,
2009
, “
Viscoelastic Relaxation and Recovery of Tendon
,”
Ann. Biomed. Eng.
,
37
(
6
), pp.
1131
1140
.10.1007/s10439-009-9687-0
24.
Fung
,
Y. C.
,
1981
,
Biomechanics: Mechanical Properties of Living Tissues
,
Springer
,
New York
.
25.
De Pascalis
,
R.
,
Abrahams
,
I. D.
, and
Parnell
,
W. J.
,
2014
, “
On Non-Linear Viscoelastic Deformations: A Reappraisal of Fung's Quasi-Linear Viscoelastic Model
,”
Proc. R. Soc. A
,
470
(
2166
).10.1098/rspa.2014.0058
26.
Balbi
,
V.
,
Shearer
,
T.
, and
Parnell
,
W. J.
,
2018
, “
A Modified Forumlation of Quasi-Linear Viscoelasticity for Transversely Isotropic Materials Under Finite Deformation
,”
Proc. R. Soc. A
,
474
(
2217
).10.1098/rspa.2018.0231
27.
De Pascalis
,
R.
,
Parnell
,
W. J.
,
Abrahams
,
I. D.
,
Shearer
,
T.
,
Daly
,
D. M.
, and
Grundy
,
D.
,
2018
, “
The Inflation of Viscoelastic Balloons and Hollow Viscera
,”
Proc. R. Soc. A
,
474
(
2218
).10.1098/rspa.2018.0102
28.
Pipkin
,
A. C.
, and
Rogers
,
T. G.
,
1968
, “
A Non-Linear Integral Representation for Viscoelastic Behaviour
,”
J. Mech. Phys. Solids
,
16
(
1
), pp.
59
72
.10.1016/0022-5096(68)90016-1
29.
Woo
,
S. L.-Y.
,
Gomez
,
M. A.
, and
Akeson
,
W. H.
,
1981
, “
The Time and History-Dependent Viscoelastic Properties of the Canine Medial Collateral Ligament
,”
ASME J. Biotech. Eng.
,
103
(
4
), pp.
293
298
.10.1115/1.3138295
30.
Haut
,
R. C.
, and
Little
,
R. W.
,
1972
, “
A Constitutive Equation for Collagen Fibers
,”
J. Biotech
,
5
(
5
), pp.
425
430
.10.1016/0021-9290(72)90001-2
31.
DeFrate
,
L. E.
, and
Li
,
G.
,
2007
, “
The Prediction of Stress-Relaxation of Ligaments and Tendons Using the Quasi-Linear Viscoelastic Model
,”
Biotech. Model. Mechanobiol.
,
6
(
4
), pp.
245
251
.10.1007/s10237-006-0056-8
32.
Pioletti
,
D. P.
, and
Rakotomanana
,
L. R.
,
2000
, “
On the Independence of Time and Strain Effects in the Stress Relaxation of Ligaments and Tendons
,”
J. Biotech.
,
33
(
12
), pp.
1729
1732
.10.1016/S0021-9290(00)00128-7
33.
Limbert
,
G.
, and
Middleton
,
J.
,
2004
, “
A Transversely Isotropic Viscohyperelastic Material: Application to the Modeling of Biological Soft Connective Tissues
,”
Int. J. Solids Struct.
,
41
(
15
), pp.
4237
4260
.10.1016/j.ijsolstr.2004.02.057
34.
Limbert
,
G.
, and
Middleton
,
J.
, in
2005
,
IUTAM Symposium on Impact Biomechanics: From Fundamental Insights to Applications
, An Anisotropic Viscous Hyperelastic Constitutive Model of the Posterior Cruciate Ligament Suitable for High Loading Rate Situations: Theoretical Framework and Material Identification,
M. D.
Gilchrist
, ed.,
Springer
,
Dortrecht, The Netherlands
, pp.
477
483
.
35.
Johnson
,
G. A.
,
Livesay
,
G. A.
,
Woo
,
S. L.-Y.
, and
Rajagopal
,
K. R.
,
1996
, “
A Single Integral Finite Strain Viscoelastic Model of Ligaments and Tendons
,”
ASME J. Biotech. Eng.
,
118
(
2
), pp.
221
226
.10.1115/1.2795963
36.
Decreamer
,
W. F.
,
Maes
,
M. A.
,
Vanhuyse
,
V. J.
, and
Vanpeperstraete
,
P.
,
1980
, “
A Non-Linear Viscoelastic Constitutive Equation for Soft Biological Tissues, Based Upon a Structural Model
,”
J. Biotech.
,
13
, pp.
559
564
.10.1016/0021-9290(80)90056-1
37.
Ciarletta
,
P.
, and
Ben Amar
,
M.
,
2009
, “
A Finite Dissipative Theory of Temporary Interfibrillar Bridges in the Extracellular Matrix of Ligaments and Tendons
,”
J. Roy. Soc. Interface
,
9
, pp.
909
924
.10.1098/rsif.2008.0487
38.
Thornton
,
G. M.
,
Oliynyk
,
A.
,
Frank
,
C. B.
, and
Shrive
,
N. G.
,
1997
, “
Ligament Creep Cannot Be Predicted From Stress Relaxation at Low Stress: A Biomechanical Study of the Rabbit Medial Collateral Ligament
,”
J. Orthop. Res.
,
15
(
5
), pp.
652
656
.10.1002/jor.1100150504
39.
Thornton
,
G. M.
,
Frank
,
C. B.
, and
Shrive
,
N. G.
,
2001
, “
Ligament Creep Behavior Can Be Predicted From Stress Relaxation by Incorporating Fiber Recruitment
,”
J. Rheol.
,
45
(
2
), pp.
493
507
.10.1122/1.1343877
40.
De Vita
,
R.
, and
Slaughter
,
W. S.
,
2006
, “
A Structural Constitutive Model for the Strain Rate-Dependent Behavior of Anterior Cruciate Ligaments
,”
Int. J. Solids Struct.
,
43
(
6
), pp.
1561
1570
.10.1016/j.ijsolstr.2005.04.022
41.
Sopakayang
,
R.
, and
Vita
,
R. D.
,
2011
, “
A Mathematical Model for Creep, Relaxation and Strain Stiffening in Parallel-Fibered Collagenous Tissues
,”
Med. Eng. Phys.
,
33
(
9
), pp.
1056
1063
.10.1016/j.medengphy.2011.04.012
42.
Reese
,
S. P.
, and
Weiss
,
J. A.
,
2014
, “
Tendons and Ligaments: Current State and Future Directions
,”
Multiscale Modeling in Biomechanics and Mechanobiology
,
S.
De
,
W.
Hwang
, and
E.
Kuhl
, eds.,
Springer-Verlag
,
London, UK
, pp.
159
206
43.
Thompson
,
M. S.
,
Bajuri
,
M. N.
,
Khayyeri
,
H.
, and
Isaksson
,
H.
,
2017
, “
Mechanbiological Modelling of Tendons: Review and Future Opportunities
,”
Proc. Inst. Mech. E, Part H
,
231
(
5
), pp.
369
377
.10.1177/0954411917692010
44.
Lanir
,
Y.
,
1980
, “
A Microstructure Model for the Rheology of Mammalian Tendon
,”
ASME J. Biotech. Eng.
,
102
(
4
), pp.
332
339
.10.1115/1.3138231
45.
Sverdlik
,
A.
, and
Lanir
,
Y.
,
2002
, “
Time-Dependent Mechanical Behaviour of Sheep Digital Tendons, Including the Effects of Pre-Conditioning
,”
ASME J. Biotech. Eng.
,
124
(
1
), pp.
79
84
.10.1115/1.1427699
46.
Raz
,
E.
, and
Lanir
,
Y.
,
2009
, “
Recruitment Viscoelasticity of the Tendon
,”
ASME J. Biotech. Eng.
,
131
(
11
), p.
111008
.10.1115/1.3212107
47.
Gupta
,
H. S.
,
Seto
,
J.
,
Krauss
,
S.
,
Boesecke
,
P.
, and
Screen
,
H. R. C.
,
2010
, “
In Situ Multi-Level Analysis of Viscoelastic Deformation Mechanisms in Tendon Collagen
,”
J. Struct. Biol.
,
169
(
2
), pp.
183
191
.10.1016/j.jsb.2009.10.002
48.
Screen
,
H. R. C.
,
Seto
,
J.
,
Krauss
,
S.
,
Boesecke
,
P.
, and
Gupta
,
H. S.
,
2011
, “
Extrafibrillar Diffusion and Intrafibrillar Swelling at the Nanoscale Are Associated With Stress Relaxation in the Soft Collagenous Matrix Tissue of Tendons
,”
Soft Matter
,
7
(
23
), pp.
11243
11251
.10.1039/c1sm05656e
49.
Puxkandl
,
R.
,
Zizak
,
I.
,
Paris
,
O.
,
Keckes
,
J.
,
Tesch
,
W.
,
Bernstorff
,
S.
,
Purslow
,
P.
, and
Fratzl
,
P.
,
2002
, “
Viscoelastic Properties of Collagen: Synchrotron Radiation Investigations and Structural Model
,”
Philos. Trans. R. Soc. London B
,
357
(
1418
), pp.
191
197
.10.1098/rstb.2001.1033
50.
Screen
,
H. R. C.
,
Lee
,
D. A.
,
,
D. L.
, and
Shelton
,
J. C.
,
2004
, “
An Investigation Into the Effects of the Hierarchical Structure of Tendon Fascicles on Micromechanical Properties
,”
J. Eng. Med.
,
218
(
2
), pp.
109
119
.10.1243/095441104322984004
51.
Legerlotz
,
K.
,
Jones
,
G. C.
,
Screen
,
H. R. C.
, and
Riley
,
G. P.
,
2013
, “
,”
Scand. J. Med. Sci. Sport
,
23
(
1
), pp.
31
37
.10.1111/j.1600-0838.2011.01410.x
52.
Baldwin
,
S. J.
,
Quigley
,
A. S.
,
Clegg
,
C.
, and
Kreplak
,
L.
,
2014
, “
Nanomechanical Mapping of Hydrated Rat Tail Tendon Collagen I Fibrils
,”
Biophys. J.
,
107
(
8
), pp.
1794
1801
.10.1016/j.bpj.2014.09.003
53.
Wenger
,
M. P. E.
,
Bozec
,
L.
,
Horton
,
M. A.
, and
Mesquida
,
P.
,
2007
, “
Mechanical Properties of Collagen Fibrils
,”
Biophys. J.
,
93
(
4
), pp.
1255
1263
.10.1529/biophysj.106.103192
54.
Grant
,
C. A.
,
Phillips
,
M. A.
, and
Thomson
,
N. H.
,
2012
, “
Dynamic Mechanical Analysis of Collagen Fibrils at the Nanoscale
,”
J. Mech. Behav. Biomed. Mater.
,
5
(
1
), pp.
165
170
.10.1016/j.jmbbm.2011.08.020
55.
Henninger
,
H. B.
,
Valdez
,
W. R.
,
Scott
,
S. A.
, and
Weiss
,
J. A.
,
2015
, “
Elastin Governs the Mechanical Response of Medial Collateral Ligament Under Shear and Transverse Tensile Loading
,”
Acta Biomater.
,
25
, pp.
304
312
.10.1016/j.actbio.2015.07.011
56.
Shearer
,
T.
,
Thorpe
,
C. T.
, and
Screen
,
H. R. C.
,
2017
, “
The Relative Compliance of Energy-Storing Tendons May Be Due to the Helical Fibril Arrangement of Their Fascicles
,”
J. R. Soc. Interface
,
14
(
133
).10.1098/rsif.2017.0261
57.
Screen
,
H. R. C.
,
Shelton
,
J. C.
,
Chhaya
,
V. H.
,
Kayser
,
M. V.
,
,
D. L.
, and
Lee
,
D. A.
,
2005
, “
The Influence of Noncollagenous Matrix Components on the Micromechanical Environment of Tendon Fascicles
,”
Ann. Biomed. Eng.
,
33
(
8
), pp.
1090
1099
.10.1007/s10439-005-5777-9
58.
Thorpe
,
C. T.
,
Udeze
,
C. P.
,
Birch
,
H. L.
,
Clegg
,
P. D.
, and
Screen
,
H. R. C.
,
2012
, “
Specialization of Tendon Mechanical Properties Results From Interfascicular Differences
,”
J. R. Soc. Interface
,
9
(
76
), pp.
3108
3117
.10.1098/rsif.2012.0362
59.
Bevan
,
T.
,
Merabet
,
N.
,
Hornsby
,
J.
,
Watton
,
P. N.
, and
Thompson
,
M. S.
,
2018
, “
A Biomechanical Model for Fibril Recruitment: Evaluation in Tendons and Arteries
,”
J. Biotech.
,
74
, pp.
192
196
.10.1016/j.jbiomech.2018.03.047
60.
Khayyeri
,
H.
,
Gustafsson
,
A.
,
Heuijerjans
,
A.
,
Matikainen
,
M. K.
,
Julkunen
,
P.
,
Eliasson
,
P.
,
Aspenberg
,
P.
, and
Isaksson
,
H.
,
2015
, “
A Fibre-Reinforced Poroviscoelastic Model Accurately Describes the Biomechanical Behaviour of the Rat Achilles Tendon
,”
PLoS One
,
10
(
6
), p.
e0126869
.10.1371/journal.pone.0126869
61.
Roy
,
S.
,
Boss
,
C.
,
Rezakhaniha
,
R.
, and
Stergiopulos
,
N.
,
2010
, “
Experimental Characterization of the Distribution of Collagen Fiber Recruitment
,”
J. Biorheol.
,
24
(
2
), pp.
84
93
.10.1007/s12573-011-0027-2
62.
Shearer
,
T.
,
Rawson
,
S.
,
Castro
,
S. J.
,
Ballint
,
R.
,
,
R. S.
,
Lowe
,
T.
,
Vila-Comamala
,
J.
,
Lee
,
P. D.
, and
Cartmell
,
S. H.
,
2019
, “
X-Ray Computed Tomography of the Anterior Cruciate Ligament and Patellar Tendon
,”
Muscles Ligaments Tendons J.
,
04
(
02
), pp.
238
244
.10.32098/mltj.02.2014.26
63.
Balint
,
R.
,
Lowe
,
T.
, and
Shearer
,
T.
,
2016
, “
Optimal Contrast Agent Staining of Ligaments and Tendons for X-Ray Computed Tomography
,”
PLoS One
,
11
(
4
), p.
e0153552
.10.1371/journal.pone.0153552
64.
Shearer
,
T.
,
,
R. S.
,
Hidalgo-Bastida
,
A.
,
Sherratt
,
M. J.
, and
Cartmell
,
S. H.
,
2016
, “
Three-Dimensional Visualisation of Soft Biological Structures by X-Ray Computer Micro-Tomography
,”
J. Cell Sci.
,
129
(
13
), pp.
2483
2492
.10.1242/jcs.179077