A new strain energy function for the hyperelastic modelling of ligaments and tendons based on fascicle microstructure

A new strain energy function for the hyperelastic modelling of ligaments and tendons based on the geometrical arrangement of their fibrils is derived. The distribution of the crimp angles of the fibrils is used to determine the stress-strain response of a single fascicle, and this stress-strain response is used to determine the form of the strain energy function, the parameters of which can all potentially be directly measured via experiments - unlike those of commonly used strain energy functions such as the Holzapfel-Gasser-Ogden (HGO) model, whose parameters are phenomenological. We compare the new model with the HGO model and show that the new model gives a better match to existing stress-strain data for human patellar tendon than the HGO model, with the average relative error in matching this data when using the new model being 0.053 (compared with 0.57 when using the HGO model), and the average absolute error when using the new model being 0.12 MPa (compared with 0.31 MPa when using the HGO model).


Introduction
Ligaments and tendons are fundamental structures in the musculoskeletal systems of vertebrates. Ligaments connect bone to bone to provide stability and allow joints to function correctly, whereas tendons connect bone to muscle to allow the transfer of forces generated by muscles to the skeleton. The wide variety of roles played by different ligaments and tendons requires them to have considerably different mechanical responses to applied forces, and their differing stress-strain behaviours have been well documented (Benedict et al., 1968;Tipton et al., 1986).
Ligaments and tendons consist of collagenous fibres organised in a hierarchical structure (Kastelic et al., 1978). Their main subunit is the fascicle, which consists of fibrils arranged in a crimped pattern (see Fig. 1). Further subunits in the hierarchy can be observed; however, the mechanics on these lengthscales will not be considered here. Instead, we shall focus on the effect of the geometrical arrangement of the fibrils within fascicles on the stress-strain properties of ligaments and tendons.
From a modelling perspective, ligaments and tendons can be categorised as fibre-reinforced biological soft tissues. A wide variety of models has been proposed to describe such tissues; however, to the author's knowledge, there has not yet been a successful attempt to develop a constitutive model within a non-linear elastic framework that includes the required anisotropy and characteristic stress-strain behaviour, which is non-linear with increasing stiffness for small strains (this region of the stress-strain curve is commonly termed the toe region) and subsequently linear, and, crucially, depends only on directly measurable parameters. Existing models are either phenomenological (Fung, 1967;Holzapfel et al., 2000), or lacking in the required material properties (such as the neo-Hookean model, which was developed for modelling rubber, but has still been used extensively in modelling biological soft tissues Miller, 2001Miller, , 2005, or both (Gou, 1970).
Early work on modelling biological tissue was carried out by Fung (1967). Fung showed that the stress in rabbit mesentery under uniaxial tension appears to increase exponentially as a function of increasing stretch. This exponential stress-strain relationship appears to approximate the behaviour of many biological soft tissues well, but only in a phenomenological sense and there is no microstructural basis for the choice of the exponential function. In 1970, Gou built upon Fung's work and proposed an isotropic strain energy function (SEF) for biological tissues that similarly gives an exponential stressstrain relationship in the case of uniaxial tension, but since this model is isotropic, it is not suitable for modelling anisotropic tissues such as ligaments and tendons.
With regard specifically to ligaments and tendons, various models were proposed over the following decades, as summarised in the review article by Woo et al. (1993). The models proposed involved infinitesimal elasticity (Frisen et al., 1969), finite elasticity (Hildebrant Contents (Fung, 1968) and single integral finite strain viscoelasticity theory (Johnson et al., 1992). In particular, we note the work of Kastelic et al. (1980), in which a model was developed for the stress-strain response of a fascicle, taking into account the distribution of fibril crimp. It was shown that a radial variation in the crimp angle of a fascicle's fibrils could lead to a nonlinear stress-strain relationship of the form typically observed in tension tests. Unfortunately, however, an error in the implementation of Hooke's law in that paper led to the derived relationship being incorrect, as we discuss further in Section 2.
Arguably the most influential model to be developed in the last 20 years for modelling biological tissues is the SEF proposed by Holzapfel et al. (2000), often referred to as the Holzapfel-Gasser-Ogden (HGO) model: where I 1 and I 4 are strain invariants, defined by where C ¼ F T F is the right Cauchy-Green tensor, where F is the deformation gradient tensor (Ogden, 1997), and M is a unit vector pointing in the direction of the tissue's fibres before any deformation has taken place, c, k 1 and k 2 are material parameters, and the above expression is only valid when I 4 Z 1 (when I 4 o 1, W ¼ c=2ðI 1 À 3Þ). This SEF was proposed as a constitutive model for arteries and is commonly used, along with its variants (Holzapfel and Ogden, 2010) to model a wide variety of biological materials. The advantages of this model are clearit retains an elegant mathematical simplicity, whilst also providing the required anisotropy and "exponentialshaped" stress-strain curve common to many biological materials; however, as it is a phenomenological model, the parameters c, k 1 and k 2 cannot be directly linked to measurable quantities, and therefore the model has restricted predictive capabilities. A large number of phenomenological, transversely isotropic, nonlinear elastic models of biological soft tissues have been proposed. The following models were collated by Murphy (2013), where the parameters c i , i¼1,2,3,4,5,6,7 are material parameters that can be chosen to match experimental data. Humphrey and Lin (1987) proposed this strain energy function for modelling passive cardiac tissue: Humphrey et al. (1990) proposed the following for passive myocardium: 1=2 4 À 1Þ 3 þ c 3 ðI 1 À 3Þþc 4 ðI 1=2 4 À 1ÞðI 1 À 3Þþc 5 ðI 1 À 3Þ 2 : Fung et al. (1993) proposed  Kastelic et al., 1978).
where T ¼ c 2 ðI 1 À 3Þ 2 þ c 3 ðI 4 À 1Þ 2 þ c 4 ðI 1 À 3ÞðI 4 À 1Þ; to model canine thoracic aorta tissue, and Chui et al. (2007) used a similar function to model porcine liver tissue: W ¼ c 1 log ð1 À TÞþc 5 ðI 1 À 3Þ 2 þ c 6 ðI 4 À 1Þ 2 þc 7 ðI 1 À 3ÞðI 4 À 1Þ: Taber (2004) used to model cardiac papillary muscles and the embryonic heart, and the simpler function has been used by many authors, for example, Destrade et al. (2008), Ning et al. (2006), and Rohrle and Pullan (2007). Weiss et al. (1996) proposed a strain energy function whose anisotropic component was expressed in terms of derivatives with respect to the fascicle stretch λ: where W λ represents the derivative of W with respect to λ and λ n is the critical stretch at which the non-linear region of the stressstrain curve ends. This strain energy function was then utilised by Peña et al. (2006) to investigate the effect of initial strains on the properties of biological tissues. All of the above strain energy functions have been shown to agree with experimental data reasonably well. The issue common to all of them, however, is that the material parameters c i , i¼ 1,2,3,4,5,6,7 have no direct physical interpretation, and therefore cannot be measured directly.
In this paper, we develop a model specifically for ligaments and tendons, but which will potentially be adaptable to other fibrous soft tissues. All parameters in this model can be directly measured via experiments. We neglect viscoelastic effects such as strain-rate dependence and hysteresis and therefore expect our model to be valid only in the low strain-rate regime where hysteretic effects are minimised. The strain energy function derived here, however, can be used to model any elastic deformation, and could potentially be incorporated into finite-strain viscoelastic models in the future. In its current form, it can easily be utilised in any finite element software that implements transversely isotropic hyperelasticity simply by replacing the existing strain energy function with that derived here.
In Section 2, we follow the method outlined by Kastelic et al. (1980) to derive stress-strain relationships for fascicles with differing fibril crimp angle distributions. By correcting the aforementioned mistake in that paper and assuming that the crimp angle distribution has a different functional form to that utilised by Kasetlic et al., we derive a simple stress-strain relationship that is used to determine the required form of the new SEF in Section 3. In Section 4, the ability of the new SEF to match experimental data is compared with that of the HGO model, and we conclude in Section 5.
2. The effect of fibril crimp angle distribution Kastelic et al. (1980) derived the stress-strain response of a fascicle under uniaxial tension as a function of the microstructural arrangement of its fibrils. Here, we summarise the procedure they used, and in the following section derive an SEF for a material whose fibres have the same stress-strain response as the fascicles modelled here. We want our SEF to have the simplest possible form that still incorporates microstructural information about the fascicle, so we make some simplifications to the method described by Kastelic et al. along the way. Kastelic et al. (1978) observed that fibril crimp angle varies throughout the radius of a fascicle, with the minimum crimp angle occurring at the centre of the fascicle, and the maximum crimp angle occurring at the fascicle's edge (see Fig. 2). If we assume that only fully extended fibrils contribute to the resistance of a fascicle to an applied load, then the resistance will increase as more fibrils become fully extended. The fibrils at the centre of the fascicle become taut first, then as the fascicle is extended, more and more fascicles within a circle of increasing radius contribute to its stiffness until the outer fibrils are finally taut. At this point, the "toe region" of the stressstrain curve ends, and the linear region begins.
Let us consider a fascicle with radiusã and define a nondimensional radial variable ρ ¼ρ=ã, whereρ is the (dimensional) radial variable, so that 0 r ρr1. Kastelic et al. (1980) assumed a crimp angle distribution of the form: where θ p ðρÞ is the crimp angle of the fibrils at radius ρ inside the fascicle, θ o is the crimp angle of the outermost fibrils, α and p are angle distribution parameters, and we note that θ o α ¼ θ i is the crimp angle of the fibril at the centre of the fascicle (ρ¼0).
Assuming that the stiffness of any extra-fibrillar matrix within the fascicle is negligible, the fascicle will be stress-free until at least one of its fibrils becomes taut. We will take this configuration as our reference configuration and hence choose θ so that the fascicle appears as in Fig. 3.
The strain in the fascicle direction, as a fibril at radius ρ becomes fully taut, is given by and the strain required to just straighten the outer fibrils is Note that Kastelic et al. (1980) included a "crimp blunting factor", b, in their model, which resulted in the above expressions for ϵ p ðρÞ and ϵ n being multiplied by ð1 À bÞ; however, since we are interested in deriving the simplest possible stress-strain relationship, we shall neglect this parameter.
At ϵ n , the non-linear "toe region" ends, and the linear elastic region begins. For a strain ϵ satisfying 0 r ϵrϵ n , however, there is, according to the model, an internal circular area of taut fibrils, each carrying a share of the applied load. The radius R p of this circle can be determined from (13), by equating ϵ ¼ ϵ p ðR p Þ: Upon using (12) in (15), and solving for R p , we obtain Fibrils outside this radius retain their crimping and experience no load. As the fascicle is stretched, this radius increases until R p ¼1, at which point all fibrils are finally taut. The tensile load experienced by the fascicle P p can be determined by integration: where σ p ðρÞ is the stress in the fibrils at radius ρ, and the upper limit of integration is determined by (16).
We will assume that the fibrils obey Hooke's law (as observed by Sasaki and Odajima, 1996). We note at this point that Kastelic et al. (1980) stated Hooke's law as where the so-called "elastic deformation" Δϵ p ðρÞ was given by Δϵ p ðρÞ ¼ ϵÀϵ p ðρÞ ¼ ϵÀ This "elastic-deformation" is not the fibril strain, and as we show below, differs from the fibril strain by a quantity that is dependent on ρ.
where here ϵ f p ðρÞ is the strain in the fibrils at radius ρ, and now we can identify E as the fibril Young's modulus. The correct expression for this strain is This can be seen by considering a fibril at radius ρ of initial length l p ðρÞ within a fascicle of initial length L (see Fig. 4). We note that l p ðρÞ and L are related by l p ðρÞ cos ðθ p ðρÞÞ ¼ L: We assume that the fascicle is stretched to a point beyond which the fibril has become taut. At this point the fascicle's length is L þ ΔL, the fibril's length is l p ðρÞþΔl p ðρÞ, and We note that the strain in the fascicle is Therefore, using Eq. (23), we can write Substituting (20) and (21) into (17), we can derive an expression for the (non-dimensionalised) average traction in the direction of the fascicle: For certain values of p, (29) can be evaluated explicitly, and upon doing so for p¼ 1,2, we obtain In order to derive a strain energy function whose fibres have the same properties as a fascicle with the microstructure described here, we would be required to integrate the above expressions. The presence of the cos À 1 ðϵ þ1Þ terms in the above expressions  would make the resulting strain energy function unnecessarily complex; however, by choosing a different form for θ p ðρÞ, we can obtain a simpler expression. If, instead of (12), we choosê θ p ðρÞ ¼ sin À 1 ð sin ðθ o Þρ p Þ; whereÁ notation is used to differentiate between the new distribution and that of Kastelic et al.,then,instead of (30) and (31), we obtain To proceed, the variation in crimp angle should be measured experimentally and plotted as a function of radius, and the parameters p and θ o should be chosen such that θ p ðρÞ, orθ p ðρÞ, gives a good approximation to the observed distribution. The simple stress-strain relationship that results from usingθ p ðρÞ when p¼ 1 makes this choice of distribution function preferable; however, quantitatively, there is little difference between θ p ðρÞ andθ p ðρÞ, as can be seen in Fig. 5 where we plot θ p ðρÞ andθ p ðρÞ for p ¼ 1; 2; 3, with θ o ¼ π=6.
The maximum difference between the distributions is less than 0.4% of the value of θ p ðρÞ for this choice of θ o , and less than 1.8% of the value of θ p ðρÞ for any θ o in the range 0 r θ o r π=4.
Eqs. (30)-(34) give the contribution of the (non-dimensionalised) average traction, acting on a cross-section normal to the fascicle direction, in the direction of the fascicle, and are plotted in Fig. 6. All four expressions have the required "toe region" behaviour, and again, we see that using the new fibril distribution has a very small effect when compared with using the distribution utilised by Kastelic et al. Upon taking Taylor series expansions of these stress-strain relationships about ϵ¼0, we observe that there is no linear term. A material with the microstructure described above, therefore, should not be modelled as linear elastic near zero strain.
Clearly,τ 1 =E is the simplest of the stress-strain relationships above, so we use this expression to derive our SEF in the following section. We recall that (33) only holds for 0 r ϵrϵ n , and for ϵ4ϵ n , we havê By using the relationship λ ¼ ϵþ1, where λ is the stretch in the direction of the fascicle, we can writê This form will be used in the derivation of our SEF in the following section.

Derivation of the strain energy function
In this section, we model the ligament or tendon under consideration as an incompressible, anisotropic, hyperelastic material. For a detailed account of the relevant theory the reader is referred to Holzapfel and Ogden (2010).
We characterise our material via a SEF W, and make the commonly used assumption that W is a function of the strain invariants I 1 and I 4 , only: where ϕ is the volume fraction of the fascicles. In the current context of ligaments and tendons, W f is the strain energy associated with the fascicles and W m is the strain energy associated with the matrix they run through. We note that I 4 , which is defined in Eq. (2), can be interpreted as the square of the stretch in the fibre direction, and can be written as where m ¼ FM is the push forward of M to the deformed configuration. Given the SEF above, the Cauchy stress tensor takes the following form (this can be seen by taking ∂W=∂I 2 ¼ ∂W=∂I 5 ¼ 0 in Eq. (2.6) of Holzapfel and Ogden, 2010): where Q is a Lagrange multiplier associated with the incompressibility constraint, I is the identity tensor, and B ¼ FF T is the left Cauchy-Green tensor. In the above, T f ¼ 2W 4 m m=ϕ is the component of the Cauchy stress associated with the fascicles, which will be used to derive an expression equivalent toτ 1 in (36) and (37). The traction associated with T f acting on a face normal to the direction of the fascicles in the deformed configuration (i.e. a face with unit normal given bym ¼ m=jmj) is By equating (42) with (36) and (37), we therefore obtain two equations for the required form of W f : 2I 4 dW f dI 4 Since I 4 is the square of the stretch in the fibre direction, λ ¼ ffiffiffiffi 2I 4 dW f dI 4 ¼ Eðβ ffiffiffiffi I 4 p À 1Þ; I 4 4 1 and, hence Integrating and (47) and (48) with respect to I 4 , we obtain the required form for the anisotropic component of our SEF: where γ and η are constants of integration, which must be chosen to satisfy W f j I 4 ¼ 1 ¼ 0 and to ensure that W f is continuous at Upon applying these conditions, we obtain Finally, we note that for I 4 o 1, W f ¼0. Therefore, we now have an explicit expression for the anisotropic part of our SEF.
To determine the form of the isotropic component of our SEF, we follow the convention used by Holzapfel et al. (2000). In the modelling of arteries, a neo-Hookean component of the SEF is commonly used to model the arterial elastin which makes up the matrix and this approach is backed up by the work of Gundiah et al. (2007). We make the assumption that the loose connective tissue that forms the matrix in ligaments and tendons has similar mechanical properties to arterial elastin so that the neo-Hookean model is still reasonable. We, therefore, have where μ40 is a stress-like parameter, which for a neo-Hookean material in isolation may be identified as the ground state shear modulus.
We have now determined the complete form of our SEF, which is written explicitly as follows: where β is defined in (35), and η is defined in (51). We note that all of the parameters in the above model can be measured. The outer fibril crimp angle of a fascicle θ o can be measured via scanning electron microscopy (Yahia and Drouin, 1989), the fibril Young's modulus E can be measured via X-ray diffraction (Sasaki and Odajima, 1996), or atomic force microscopy (Svensson et al., 2012), and it should soon be possible to measure the fibre volume fraction ϕ and fibre alignment vector M via X-ray computed tomography (Shearer et al., 2014). The only parameter which is not straightforward to determine is μ; however, the techniques used by Gundiah et al. (2007) to determine the value of this parameter for arterial elastin could potentially be adapted to do the same for ligament or tendon matrix. Finally, we note that the SEF above behaves isotropically in the small strain limit. This may seem unusual for a transversely isotropic SEF, but this behaviour arises due to the microstructure of the fascicles, which do not have a linear term in their stressstrain relationship for small strains, as discussed in Section 2.

Comparison with the HGO model
In this section, we compare the new model (53)-(55) with the HGO model (1). The two models are plotted against stress-strain data for human patellar tendon taken from Johnson et al. (1994). Johnson et al. plotted separate results for samples taken from younger (aged 29-50) and older (aged 64-93) donors. We fit the models to the experiments performed on the younger samples. The patellar tendon was chosen due to the fact that its fascicles are very strongly aligned with its longitudinal axis (Shearer et al., 2014), so that the fibre alignment vector M could be chosen to coincide with the tendon's axis in the model. In this simple example, we will model the patellar tendon as a circular cylinder. This is not the case in reality, but for the purpose of illustrating the model developed above, this geometry will be sufficient.
Consider a cylinder with radius A and length B in the undeformed configuration. We describe its geometry in this configuration in terms of the polar coordinates ðR; Θ; ZÞ, so that 0 r Rr A, 0 r Θr2π, 0r Z rB. After the application of a homogeneous longitudinal stretch, its geometry in the deformed configuration is described by the polar coordinates ðr; θ; zÞ, so we have 0r r r a, 0 r θ r 2π, 0r z rb, where a and b are the deformed counterparts of A and B, respectively. The deformation, then, is described by where ζ¼b/B is the uniform stretch in the longitudinal direction and the first equation is a consequence of the fact that we are using incompressible constitutive models. Due to the symmetry of the imposed deformation, the deformation gradient tensor will be diagonal and its entries will be the principal stretches: where e i , i ¼ ðr; θ; zÞ, and E J , J ¼ ðR; Θ; ZÞ, are the deformed and undeformed unit vectors in the radial, azimuthal and longitudinal directions, respectively, and the left Cauchy-Green tensor is given by Since the fascicles of the patellar tendon are aligned with its longitudinal axis, we have M ¼ E Z , and, therefore, m ¼ ζe z . The strain invariants I 1 and I 4 can also be calculated: Upon using Eqs. (57)-(59) in Eq. (40), solving the static equilibrium equations, div T ¼ 0, and applying traction-free boundary conditions on r ¼a, we obtain explicit expressions for the longitudinal stresses T HGO zz and T zz : cðζ 2 À ζ À 1 Þþ4k 1 ζ 2 ðζ 2 À 1Þe k 2 ðζ 2 À 1Þ 2 ; ζ Z 1; where the superscript HGO is to differentiate the HGO model from that presented above. In Fig. 7, we plot the nominal stresses S HGO zz ¼ T HGO zz =ζ and S zz ¼ T zz =ζ (which give the force per unit undeformed area, to match the experimental data) as a function of the engineering strain e¼ζÀ1, along with the patellar tendon stress-strain data taken from Fig. 4 of Johnson et al. (1994). We note that this experimental data clearly displays non-linear behaviour, even for strains of less than 1%.
Since the stiffness of ligament and tendon matrix is insignificant compared with that of its fascicles, c and ð1 À ϕÞμ were chosen to be small ðc ¼ where S exp zz is the experimental stress, and the absolute error, defined at a given strain by Δ ¼ jS exp zz À S zz j; Δ HGO ¼ jS exp zz ÀS HGO zz j: The average relative error for the new model was δ ¼ 0:053, and that of the HGO model was δ HGO ¼ 0:57, and the average absolute error for the new model was Δ ¼ 0:12 MPa, and that of the HGO model was Δ HGO ¼ 0:31 MPa. Hence, the relative error associated with the new model is less than 10% of that of the HGO model and the absolute error associated with the new model is less than 41% of that of the HGO model.

Discussion
In this paper, we derive a new SEF for the hyperelastic modelling of ligaments and tendons. The form of this SEF is determined from the stress-strain response of fascicles as a function of their microstructure, using the method described by Kastelic et al. (1980). In Section 2, a simple, explicit form for this stress-strain behaviour is determined by altering the form of the fibril crimp angle distribution assumed by Kastelic et al. and the corresponding SEF is derived in Section 3. Crucially, all of the parameters in this new model can potentially be measured experimentally.
In Section 4, we compare the ability of the new model to reproduce experimental stress-strain data for human patellar tendon taken from Johnson et al. (1994) with that of the HGO model proposed in Holzapfel et al. (2000). Since independent measurements of the collagen volume fraction ϕ, fibril stiffness E, and average outer fibril crimp angle θ o are not readily available, these parameters are chosen to fit the experimental data. The new model is found to be more successful at modelling this experimental data, and the values of the parameters used in the new model are physically realistic. The value of ϕE used above gives a value of E that lies well within the range of reported values of 0.07 (Yang et al., 2008)-5.1 GPa (Cusack and Miller, 1979), assuming 0:11r ϕr1. To the author's knowledge, data on typical values of θ o is not readily available; however, the predicted value of θ o ¼ 10:71 seems reasonable. We note that the new model performs particularly well in terms of relative error. This is due to the better performance of the new model within the "toe region" of the stress-strain curve when compared with the HGO model.
In order to fully validate this model, it will be important to independently measure ϕ, E and θ o for a given ligament or tendon, input these quantities into the model and compare the predicted results with experimental stress-strain data for that ligament or tendon.

Conflict of interest statement
The author has no conflicts of interest to declare.