A new strain energy function for modelling ligaments and tendons whose fascicles have a helical arrangement of fibrils

A new strain energy function for the hyperelastic modelling of ligaments and tendons whose fascicles have a helical arrangement of fibrils is derived. The stress-strain response of a single fascicle whose fibrils exhibit varying levels of crimp throughout its radius is calculated and used to determine the form of the strain energy function. The new constitutive law is used to model uniaxial extension test data for human patellar tendon and is shown to provide an excellent fit, with the average relative error being 9.8%. It is then used to model shear and predicts that the stresses required to shear a tendon are much smaller than those required to uniaxially stretch it to the same strain level. Finally, the strain energy function is used to model ligaments and tendons whose fascicles are helical, and the relative effects of the fibril helix angle, the fascicle helix angle and the fibril crimp variable are compared. It is shown that they all have a significant effect; the fibril crimp variable governs the non-linearity of the stress-strain curve, whereas the helix angles primarily affect its stiffness. Smaller values of the helix angles lead to stiffer tendons; therefore, the model predicts that one would expect to see fewer helical sub-structures in stiff positional tendons, and more in those that are required to be more flexible.


Introduction
Ligaments and tendons are important connective tissues; ligaments connect bone to bone, providing stability and allowing joints to function correctly, and tendons connect muscle to bone to transfer forces generated by muscles to the skeleton. They both have a hierarchical structure consisting of several fibrous subunits (Kastelic et al., 1978;Screen et al., 2004), which, from largest to smallest, can be defined as follows: fascicles (50-400 μm diameter), fibrils (50-500 nm), sub-fibrils (10-20 nm), microfibrils (3-5 nm), and finally, the tropocollagen molecule (∼1.5 nm). The geometrical arrangement of many of these subunits varies between different ligaments and tendons; for example, the patellar tendon's fascicles are coaligned with its longitudinal axis, whereas the anterior cruciate ligament's are helical (Shearer et al., 2014). The fibrils within a fascicle may also either be coaligned or helical with respect to its longitudinal axis (Yahia and Drouin, 1989). In both cases, the fibrils exhibit an additional waviness, called crimp, which is superimposed upon their average direction and varies in magnitude throughout the fascicle's radius (Kastelic et al., 1978;Yahia and Drouin, 1989). This intricate structure produces complex mechanical behaviour such as anisotropy, viscoelasticity and non-linearity, which varies between different ligaments and tendons (Benedict et al., 1968;Tipton et al., 1986). It is not currently known, however, which levels of the hierarchy are most influential in governing their mechanical performance.
To begin understanding these mechanical features, it is of interest to model their elastic properties, neglecting viscoelasticity. Elastic models are expected to be valid in both the low and extremely high strain rate limits where hysteresis is minimised. Ligament and tendon stress-strain behaviour under uniaxial tension is characterised by an initial non-linear region of increasing stiffness, termed the toe-region, followed by a linear region before the onset of failure (Fig. 1). Several authors have derived expressions to describe this behaviour (Frisen et al., 1969;Kastelic et al., 1980;Kwan and Woo, 1989); however, to consider more complex deformations, it is useful to characterise the elasticity of a material in terms of a strain energy function (SEF).
Many non-linear elastic SEFs have been proposed for soft tissues (Fung, 1967;Gou, 1970;Holzapfel et al., 2000), but few have focused specifically on ligaments and tendons. Whilst many SEFs are general enough to be applied to modelling ligaments and tendons, the majority of them contain variables that cannot be directly experimentally measured (Shearer, 2015). This limits their ability to analyse which physical quantities are most important in governing a specific ligament's or tendon's behaviour. Microstuctural models are better equipped to facilitate this analysis, provided their parameters can be experimentally determined. Two examples of microstructural SEFs are those derived by Grytz and Meschke (2009) and Shearer (2015). Both are based on the geometrical arrangement of the fibrils within the fascicle and neglect any subunits below the fibril level (Fig. 2). Shearer considered a fascicle whose fibrils are coaligned with its axis, but have a distribution of crimp levels throughout its radius ( Fig. 3 (1) and (2)), whereas Grytz and Meschke considered a helical arrangement of fibrils, but neglected their crimp (Fig. 3(3) and (4)). Grytz and Meschke defined the angle that these fibrils make with the fascicle's longitudinal axis as the crimp angle, but this is not the usual definition of crimp. Here, this quantity is referred to as the fibril helix angle. The logical extension of these models is to allow the fibrils to be helically arranged and crimped ( Fig. 3(5)); this is the case considered here. A scanning electron microscope (SEM) image displaying fibrils that are both helically arranged and crimped appears in Fig. 9 of Yahia and Drouin (1989).
A considerable amount of work has been dedicated to modelling other types of fibre-reinforced composite materials. Crossley et al. (2003), for example, derived analytical solutions that govern the bending and flexure of helically reinforced, anisotropic, linear elastic cylinders. This is built on a large body of literature on modelling cables (Cardou and Jolicouer, 1997) and rope (Costello, 1978(Costello, , 1997. Adkins and Rivlin (1955) discussed finite deformations of materials that are reinforced by inextensible cords, and Spencer and Soldatos (2007) considered finite deformations of fibre-reinforced elastic solids whose fibres are capable of resisting bending. Whilst these general theories are extremely valuable for certain problems, to model the behaviour of a material with a microstructure as complex as a ligament or tendon requires a more specific model.
In this paper, a new SEF that governs the behaviour of a ligament or tendon with the microstructure described above is derived. In Section 2, the stress-strain response of a single fascicle is calculated and this relationship is used to determine the form of the new SEF in Section 3. The SEF is used to model the mechanical response of human patellar tendon to uniaxial extension and shear in Section 4. In Section 5, the case of a ligament or tendon with helical fascicles is explored and the relative effects of the fibril crimp angle, fibril helix angle and fascicle helix angle are analysed. Finally, the implications of the model are discussed in Section 6.
2. The stress-strain response of a fascicle with helically aligned fibrils Kastelic et al. (1980) derived the stress-strain response of a fascicle with fibrils that are coaligned with its longitudinal axis, and Shearer (2015) adapted their method to derive analytical expressions for these relationships for different fibril crimp angle distributions. Here, this work is extended to ligaments and tendons whose fascicles have a helical arrangement of fibrils. Kastelic et al. (1978) observed that crimp angle varies throughout the radius of a fascicle with longitudinal fibrils. It was then noted by Yahia and Drouin (1989) that this is also the case in fascicles with a helical arrangement of fibrils. The minimum crimp angle occurs at the fascicle's centre, the maximum at its edge. Therefore, assuming that only fully extended fibrils contribute to  average force per unit undeformed area acting on tendon the fascicle's resistance to an applied load, its rigidity will increase as more and more fibrils become fully extended. The fibrils at the centre of the fascicle are the first to tauten, then as the fascicle is increasingly stretched, more and more fibrils within a circle of increasing radius begin to contribute to its stiffness until the outer fibrils are finally fully extended. At this point, the toe-region ends and the linear region begins.
It is assumed that the crimp angle distribution is described by Shearer (2015) sin sin , 1 where ρ is a non-dimensional radial variable scaled on the fascicle's radius so that 0 1 ρ ≤ ≤ , θ o is the crimp angle of the outer fibrils and p is a parameter that can be chosen to fit an observed crimp distribution. It is assumed that the fibrils make an angle α with the fascicle's longitudinal axis, so that we can write a unit vector pointing in the fibril direction in cylindrical coordinates as follows: In the following, it is important to distinguish between the fascicle stretch λ, the component of that stretch in the fibril direction Λ, and the stretch actually experienced by a fibril, which is only greater than one once its crimp has straightened out. If the fascicle undergoes a longitudinal stretch λ and we assume that fibrils slide freely without interacting with each other, or with the extracollagenous matrix that holds them in place, then the fibril alignment vector will be modified to In reality, the fibrils are likely to interact with each other and fibril sliding is likely to affect ligament and tendon viscoelasticity (Screen, 2008); however, to keep the model as simple as possible, it is assumed that fibril interaction is negligible in the elastic case being considered here. The component of the fascicle stretch in the fibril direction is This stretch straightens out the fibrils' crimp until they become taut. At this point, the fibrils themselves begin to stretch. The stretch in the fibril direction required to straighten the fibrils at radius ρ is 1 cos and that required to straighten the outer fibrils is Using Eq. (4), the fascicle stretch corresponding to Λ ⁎ can be calculated 1 cos , the toe-region ends. For a fascicle strain 1 λ ϵ = − satisfying 0 ≤ ϵ ≤ ϵ ⁎ , however, there is an internal area of taut fibrils, each carrying a share of the load. The radius R p of this circle can be determined from (5), by equating R Fibrils outside R p 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 is , 11 We assume that the fibrils obey a linear stress-strain relationship (as observed by Sasaki and Odajima, 1996) E , 13 where E is the fibril Young's modulus, and the fibril strain p f ρ ϵ ( ) for a given fibril within the radius R p is (Shearer, 2015) cos 1. 14 Upon substituting (13) and (14) into (12) The critical stretch λ ⁎ at which the toe-region ends is a monotonically increasing function of α and θ o , provided 0 and there is no toeregion. This makes sense since, if 0 o θ = , all the fibrils have no crimp and contribute to the tendon's stiffness from the start of its deformation.

Derivation of the strain energy function
The ligament or tendon under consideration is modelled as incompressible, anisotropic and hyperelastic and is characterised via the SEF W. It is assumed that W is a function of the strain invariants I 1 and I 4 only where ϕ is the collagen volume fraction. In the current context, W f is the strain energy associated with the fascicles and W m is that associated with the extracollagenous matrix. The invariants I 1 and I 4 are defined by is the right Cauchy-Green tensor (F is the deformation gradient, Ogden, 1997), and M is a unit vector pointing in fascicle direction in the undeformed configuration. I 4 can be interpreted as the square of the stretch in the fascicle direction where Q is a Lagrange multiplier associated with the incompressibility constraint, I is the identity tensor, is the component of the Cauchy stress associated with the fascicles, which will be used to derive an expression equivalent to τ 1 in (16) and (18).  (24) with (16) and (18) Finally, note that for I 1 4 < , W f ¼0. Therefore, the anisotropic part of the SEF is now explicit.
As in Shearer (2015), the isotropic component of the SEF is chosen to be neo-Hookean where 0 μ > is the ground state shear modulus of the extracollagenous matrix. The neo-Hookean SEF accurately models the behaviour of arterial elastin (Gundiah et al., 2007) and it is assumed that the extracollagenous matrix in ligaments in tendons has similar mechanical properties.
The complete form of the SEF can now be written explicitly This is the SEF for a ligament or tendon whose fascicles have longitudinal fibrils (Shearer, 2015).

Unixial extension
The SEF derived in the previous section is now used to determine a stress-strain relationship, which is compared with the experimental data for younger human patellar tendons from (Johnson et al., 1994). For simplicity, the tendon is modelled as a circular cylinder.
We assume that the tendon has radius A and length L in the undeformed configuration. Its geometry is described in terms of the cylindrical coordinates R Z , , After applying a homogeneous longitudinal stretch, the corresponding deformed coordinates are r z , , θ ( ), so r a 0 ≤ ≤ , 0 2 θ π ≤ < , z l 0 ≤ ≤ , where a and l are the deformed counterparts of A and L. The deformation, then, is described by is the longitudinal stretch, and the first equation is a consequence of the tendon's incompressibility. The deformation gradient is The static equations of equilibrium in the absence of body forces are where div is the divergence operator in the deformed configuration. In this case, (49) reduces to Therefore, Q is a constant, which can be determined by applying a traction-free boundary condition on r ¼ a ( T n 0 · = , where n e r = ), which gives T T T 0 rr r a r r a rz r a The first of these gives The corresponding nominal stress, which gives the force per unit undeformed area (the correct measure of stress to use when comparing with experimental data), is given by S T / zz zz ζ = . In Shearer (2015), S zz (but with α¼0°, since the patellar tendon's fascicles were assumed to have longitudinal fibrils) was fitted to the data from Johnson et al. (1994). Since the stiffness of a ligament's or tendon's matrix is negligible compared to that of its fascicles, 1 ϕ μ ( − ) was chosen to be small ( 1 0 . 0 1M P a ϕ μ ( − ) = ) relative to the reported values of E, which range from 32 MPa (Graham et al., 2004) to 11 GPa (Wenger et al., 2007). The remaining parameters were determined using the FindFit function in Mathematica 7 (Wolfram Research, Inc., Champaign, IL, 2008)  find a locally optimal solution in the non-linear case (which is the case here); therefore, the restrictions (53) are crucial to ensure that the given solution is physically realistic. Using the default settings, it attempts to make the numerical error in a result of size x be less than

Shear
Whilst ex vivo uniaxial extension data collected in the laboratory provides useful information about the mechanical behaviour of a ligament or tendon, in vivo loading conditions are more complex and are likely to involve a combination of both extension and shear. Shear tests can be extremely challenging to undertake due to the difficulties involved in gripping samples. Additionally, to obtain a cubic sample suitable for a shear test, the specimen must be dissected, which could potentially affect its mechanical properties due to damage to its ultrastructure. Mathematical modelling allows us to simulate a shear test, whilst avoiding any experimental difficulties.
Here, two distinct modes of simple shear are considered: those in which planes parallel to the fascicle direction remain undeformed:  Using the parameter values derived for human patellar tendon, above, these shear stresses are plotted in Fig. 5. As can be seen, the model predicts that: the response of the patellar tendon to shear in which planes parallel to the fascicle length remain undeformed is linear, it experiences strain stiffening under shear in which planes perpendicular to the fascicle direction remain undeformed, in both cases, the stresses required to shear the tendon are much smaller than those required to uniaxially stretch it to the same level of strain.

Ligaments and tendons with helical fascicles
In addition to having fascicles with helically arranged fibrils, some ligaments and tendons such as the anterior cruciate ligament (Shearer et al., 2014) and the extensor carpi ulnaris tendon (Kalson et al., 2011) also have fascicles that form a helix with their longitudinal axis. To model this, the fascicle alignment vector M must be modified so that it no longer only points in the Z-direction, but also has a component in the Θ-direction where ψ is the angle that the fascicles make with the ligament's or tendon's longitudinal axis. Upon doing this, and again assuming the ligament or tendon undergoes a stretch ζ, we obtain a mod-   and E 1027 MPa ϕ = , as before) is plotted as a function of the engineering strain e 1 ζ = − for several values of these parameters. In each plot, we use 20 o α ψ θ = = =°as a base case (the blue lines) and investigate the effect of changing each parameter in turn by plotting the curve when that parameter is equal to 0°(black lines), 10°(red lines), and 30°(purple lines). In  Table 1 shows the stress at 5% strain when each of the parameters is varied from the base case. It appears that θ o has the biggest effect, followed by ψ, then α; however, they all clearly have a significant effect. Note that, as mentioned in Section 2, when 0 o θ =°, there is no toe-region and the ligament or tendon behaves approximately linearly. We conclude that it is the distribution of crimp angles that primarily governs the non-linear region of a ligament's or tendon's stress-strain curve and the helix angles mainly affect their stiffness.

Discussion
In this paper, a new SEF for the hyperelastic modelling of ligaments and tendons that consist of fascicles with helically arranged fibrils has been derived. This is an extension of the work of Shearer (2015) and Grytz and Meschke (2009). Grytz and Meschke included a fibril bending stiffness in their model (i.e. they modelled the fibrils as helical rods). Whilst the inclusion of this parameter may seem sensible, in their case it led to an SEF that cannot be expressed explicitly in terms of strain invariants, and instead has to be determined algorithmically (although the inclusion of bending stiffness does not always lead to an intractable SEF, Spencer and Soldatos, 2007). By neglecting bending stiffness, we have essentially modelled the fibrils as strings, which, given the extreme aspect ratios exhibited by fibrils (Trotter and Koob, 1989), seems to be a reasonable assumption. By doing this, it has been possible to incorporate fibril crimping into the model and obtain an SEF that is explicitly expressible in terms of strain   Table 1 The effect of varying the three parameters α, ψ, θ o on the stress at 5% strain.

Parameter values
Stress ( invariants, and whose material parameters can all be directly experimentally measured. Examples of ways to measure most of these parameters are discussed in Shearer (2015). The only parameter not mentioned in that paper is the fibril helix angle, which can be measured via SEM (Yahia and Drouin, 1989). The SEF derived here is dependent on only two invariants, but at least three are needed to fully capture the mechanical response of incompressible, transversely isotropic materials under complex loading conditions (Destrade et al., 2013). A natural extension of the SEF derived here could be obtained by replacing the neo-Hookean matrix component with a Mooney-Rivlin form.
The errors in fitting the experimental data in Section 4.1 are actually slightly higher than those achieved by the original model proposed by Shearer (2015) ( 5.3% δ = , 0.12 MPa Δ = ). This may indicate that the fibril helix angle in human patellar tendon is closer to 0°than the 27°that was assumed based on images of canine patellar tendon. To fully test the model, all of its parameters should be independently measured for a specific ligament or tendon and used to predict its stress-strain behaviour. This prediction should then be compared with tension tests performed on that ligament or tendon.
In Section 4.2, it was shown that the model can be used to make qualitative predictions about how tendons respond to simple shear without having to dissect them, and therefore potentially damage their ultrastructure. Interestingly, the stress required to apply simple shear appears to be much smaller than that required for extension. This could have potential implications for tendon reconstruction surgery. If graft tendons are fixated such that shear forces are generated, large deformations could occur, which could potentially lead to rerupture. This situation could be avoided by assessing different potential graft geometries via finite element analysis. Since the model is expressed in terms of an SEF, it would be straightforward to implement it into finite element software to simulate in vivo deformations. The microstructural basis of the model makes it preferable to the phenomenonlogical finite strain models that have been implemented in commercial finite element packages to date.
To conclude, we consider the results of Section 5. All the angle parameters (α, ψ, θ o ) have a significant effect on the predicted mechanical properties. This may indicate that the fibrilar and fascicular structure of a ligament or tendon is optimised to suit its function. For example, one would expect to see fewer helical substructures in a stiff, positional tendon and more in those that are required to be more flexible. This prediction is supported by the experimental observations of Thorpe et al. (2012Thorpe et al. ( , 2013.

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