A Recruitment Model of Tendon Viscoelasticity That Incorporates Fibril Creep and Explains Strain-Dependent Relaxation

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 [3–5].

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 [37–41], 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.

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.

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.

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.

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.

where here, J₀, J₁, and J₂ 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).

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.

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 e_c 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.

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.

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, e_c, 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.

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 t_A 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.

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 t_B and t_D 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.

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

The upper limit of integration in Eq. (21) is t_C 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 = t_D 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 e_f. We now discuss the re-load curve and how we determine the time $tA′$ at which the stress attains a positive value once again.

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 ℓ₀ 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.

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 e_c 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.

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.

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.

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.

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.

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 l² norm of the residuals between the theoretical and experimental values at three points along the averaged relaxation curves, as illustrated in Fig. 7.

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.

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.

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 ℓ₀, 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.² The length distribution is plotted in nondimensional form as fibril length divided by initial fascicle length (i.e., $l(t)/L$⁠).

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.

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,² 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 [62–64]. 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.

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.

• 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).