An arterial constitutive model accounting for collagen content and cross-linking

It is apparent from the literature that the density of cross-links in collagenous tissue has a stiffening effect on the mechanical response of the tissue. This paper represents an initial attempt to characterize this effect on the elastic response, specifically in respect of arterial tissue. Two approaches are presented. First, a simple phenomenological continuum model with a cross-link-dependent stiffness is considered, and the influence of the crosslink density on the response in uniaxial tension is illustrated. In the second approach, a 3D model is developed that accounts for the relative orientation and stiffness of (two families of) collagen fibers and cross-links and their coupling using an invariant-based strainenergy function. This is also illustrated for uniaxial tension, and the influence of different cross-link arrangements and material parameters is detailed. Specialization of the model for plane strain is then used to show the effect of the cross-link orientation (relative to the fibers) and cross-link density on the shear stress versus the amount of shear deformation response. The elasticity tensor for the general (3D) case is provided with a view to subsequent finite element implementation. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Collagen is the most important structural protein in the body and is able to bear significant mechanical load within fibrous tissues ( Fratzl, 2008 ). In such tissues collagen forms a network together with cross-links which, from the solid mechanics point of view, contribute to the transmission of forces between the fibers of the network in both healthy and aged tissues ( Andriotis et al., 2018 ). Collagen is a hierarchical material ( Fratzl, 2008;Fratzl and Weinkamer, 2007 ), which is composed of a tropocollagen triple helix at the nanoscale, typically about 300 nm long. This tropocollagen is held together by intramolecular bonds. Aggregates of tropocollagen molecules connected together by cross-links form collagen fibrils which themselves group together to form collagen fibers. Enzymatic cross-links, which connect tropocollagen molecules at their ends and provide stability of the structure, contribute to the mechanical resistance of a fibril under tension. On the other hand non-enzymatic cross-links which can attach at any point along the length of a tropocollagen molecule can be JID: MPS [m3Gsc;August 16, 2019;16:12 ] We also consider a planar specialization of the model and illustrate it by application to a simple shear deformation, showing the effect of the cross-link orientation (relative to the fibers) and cross-link density on the shear stress versus amount of shear deformation. The reader who requires additional information on the subject of nonlinear elasticity and solid mechanics is referred to the monographs of Holzapfel (20 0 0) and Ogden (1997) , while a detailed explanation of the underlying constitutive theory for strongly anisotropic solids can be found in Spencer (1984) .

Model structure
Suppose the collagen fibers are embedded within an isotropic matrix with the volume fraction so that (1 − ) is the volume fraction of the matrix. Let the elastic properties of the matrix and fibers be described in terms of strain-energy functions iso and f , respectively. We consider iso to depend on the isotropic invariant I 1 = tr C , where C = F T F is the right Cauchy-Green tensor, F is the deformation gradient, and f , the energy associated with the fiber in the direction M in the reference configuration, to depend on I 4 = M · C M and also on a measure of the density of cross-links per unit length of the fiber in the direction M .
We consider the material to be incompressible ( J ≡ det F = 1 ) with the total elastic energy of this composite as where ρ can be thought of as the number of cross-links per unit length, subsequently referred to as the density of crosslinks with dimension 1/(length). As an example iso can take on a neo-Hookean form, while for f we could have k (ρ)(I 4 − 1) 2 / 2 , a standard reinforcing model. Herein k ( ρ) is the cross-link-dependent fiber stiffness. Note that the derivative k should be positive to reflect the increasing stiffness with increasing density of cross-links. Hence, from (1) , the Cauchy stress tensor σ can be calculated as where p is the Lagrange multiplier associated with the incompressibility constraint, b denotes the left Cauchy-Green tensor, I is the identity tensor, μ is the shear modulus of the neo-Hookean matrix and m = F M . The anisotropic term is only active if I 4 > 1, so under compression in the fiber direction, the matrix bears the stress. Because of the experimental data of fibrous tissue it is useful to represent f as an exponential, in this case given by where k 1 > 0 is a parameter with the dimension of stress, while k 2 > 0 is a dimensionless parameter. The Cauchy stress of the fibers is then denoted by σ f , i.e.
while the Cauchy stress for the matrix σ iso is Let us now consider a strip of tissue in the axial/circumferential plane with two families of fibers which are arranged symmetrically with respect to the axes, as indicated in Fig. 1 , where α is the angle between the axial direction and the direction of each family of fibers. The direction of the second fiber family is denoted by M with m = F M . According to Fig. 1 the matrix forms of M and M are where we have assumed that the collagen fibers have no out-of-plane component. A push forward gives where λ 1 and λ 2 are the principal stretches along the directions 1 (axial direction) and 2 (circumferential direction), respectively, while, from the incompressibility condition, λ 3 = λ −1 1 λ −1

. The invariant I 4 is then
and by symmetry we have M · C M = I 4 .
The Cauchy stress tensor σ = σ iso + σ f − p I then becomes which is diagonal with respect to the chosen axes (no shear stress), and hence its components are  We assume that the strip is under plane stress conditions with loads parallel to the circumferential and axial directions.

ARTICLE IN PRESS
Then the stress σ 33 is zero which allows for p to be eliminated from (10) and (11) to give where I 4 is given by (8) 2 , and λ 3 by the incompressibility condition. Thus, σ 11 and σ 22 are given in terms of λ 1 and λ 2 . By considering uniaxial stress with σ 22 = 0 , Eq. (14) can in principle be solved for λ 2 as a function of λ 1 , and then σ 11 is a function of λ 1 alone. We now focus on a special case, namely α = 0 . Hence, the two families of fibers coincide and the direction of the collagen fibers is the axial direction. For this case we choose λ = λ 1 , so by symmetry λ 2 = λ 3 = λ −1 / 2 , and according to (8) 2 we have I 4 = λ 2 . We use the dimensionless quantities σ 11 = σ 11 /μ and k 1 = k 1 /μ, and obtain from (13) , which is an explicit expression for σ 11 in terms of λ, where , k 1 (ρ) and k 2 need to be specified.
The functional dependence of k 1 on ρ can be modeled by any suitable function but for simplicity of illustration we consider the quadratic equation k 1 = k 0 + ā ρ 2 , where k 0 ≥ 0 is the value of k 1 at ρ = 0 and ā > 0 is a parameter with dimension (length) 2 . As ρ increases the response becomes stiffer.

A model with aligned collagen fibers connected by separately aligned cross-links
We consider collagen fibers to be in the direction of the unit vector E 1 and let E R be the radial unit vector normal to that direction. In addition we consider two symmetrically disposed families of cross-links in the directions of the unit vectors L + and L − ( L stands for link), which are defined by where α 0 defines their orientation relative to the collagen fiber direction; see  A uniaxial deformation with stretch λ is applied along the collagen fibers so that by symmetry and considering the material to be incompressible the deformation gradient has the form We define e = F E 1 = λE 1 and e r = F E R = λ −1 / 2 E R as the push-forwards of E 1 and E R under the deformation. The corresponding push-forwards of L + and L − are then For this special deformation the isotropic invariant is given by I 1 = tr C = λ 2 + 2 λ −1 , while the squares of the stretches in the directions E 1 , L + and L − are wherein the invariants I 4 and I are defined. For the cross-links to be extended, i.e. when λ > 1, α 0 has to be restricted according to which is satisfied for all λ > 1 if cos α 0 > 1 / √ 3 . The invariant I 4 is the square of the stretch in the collagen fibers and I is the square of the stretch in each of the cross-link directions. We also define the coupling between the collagen fiber and cross-link directions by the quantities I + 8 and I − 8 , which are given by These are not themselves invariants (their sign changes under reversal of either e or l ± ), but (I + 8 ) 2 = (I − 8 ) 2 is invariant. The values of the invariants I 1 , I 4 , and I and the quantities I ± 8 in the reference configuration are 3, 1, 1 and ± cos α 0 , respectively.
We now consider a strain-energy function , which is a function of I 1 , I 4 , I , I + 8 and I − 8 . Specifically we consider to have the form where and denote the volume fractions of the collagen fibers and the cross-links, respectively. The functions iso , f and c are the energies stored in the matrix material, the collagen fibers and the cross-links, respectively, while the two fc -terms represent the interaction energies between the collagen fibers and the cross-links.
Noting that fc (I − 8 ) = − fc (I + 8 ) , the Cauchy stress tensor σ can be written in the form where we have used the abbreviations Now, for the subsequent component forms we need Hence, the stress components are where the argument of fc is I + 8 . By eliminating the Lagrange multiplier p by subtraction of (27) from (26) we obtain the uniaxial stress σ = σ 11 as Now let us consider some specific energy functions. For the matrix material we use the isotropic neo-Hookean material iso = 1 2 where the constant μ is a positive parameter, and for the collagen fibers we use the standard exponential form, i.e.
where k 1 > 0 is a stress-like constant and k 2 > 0 is a dimensionless constant. There is very little if any information available about the mechanical properties of cross-links. Therefore, for simplicity of illustration, we make the assumption that c has the quadratic reinforcing form where ν is a positive parameter with the dimension of stress that measures the strength of the cross-links, and is referred to as the cross-link parameter. Note that the constant k 1 in (30) is different from that in (3) , and we now write it as k 0 temporarily. Then, by comparing the quadratic approximation of (3) with (30) and (31) we obtain ν + k 0 = k 1 (ρ) , which relates the cross-link stiffness ν to the cross-link density of the first model. Similarly to (31) , for fc we take the form fc = 1 2 where κ is also a positive stress-like parameter. It measures the strength of the interaction between the fibers and the cross-links. Hence, by using (24) and according to (28) , the Cauchy stress σ has the form +4 ν(λ 2 cos 2 α 0 + λ −1 sin 2 α 0 − 1)(λ 2 cos 2 α 0 − λ −1 sin 2 α 0 ) + 4 κ (λ 2 − 1) λ 2 cos 2 α 0 .
In Fig. 4 we plot the dimensionless stress σ = σ /μ against the stretch λ for a representative selection of the parameters involved in (33) . Fig. 4 (a) shows how the response depends on the orientation α 0 of the cross-links for a fixed value of ν = ν/μ, while Fig. 4 (b) illustrates the dependence on ν for a fixed value of α 0 , in each case for fixed values of the other parameters, as specified in the caption of Fig. 4 . It is clear that the cross-links stiffen the response. In Fig. 4 (a) the response becomes stiffer as the cross-links become more aligned with the fibers, much stiffer than in the absence of cross-links, while Fig. 4 (b) shows that an increase in the density of the cross-links likewise stiffens the response.  4. Plots of the dimensionless Cauchy stress σ = σ /μ versus the stretch λ: (a) for three values of the cross-link angle α 0 ( π/16, π/6, π/4) compared with the plot for the case of no cross-links. On the basis of (33) the following parameters were used

ARTICLE IN PRESS
for four values of the dimensionless cross-link parameter ν (10.0, 5.0, 1.0, 0). On the basis of (33) the following parameters

Formulation for a general fiber direction
In the previous section we considered a special deformation with the fiber directions aligned with the axis E 1 of extension. In the present section we generalize this for an arbitrary fiber direction and the corresponding cross-links. Consider the axis E 1 and the associated rectangular Cartesian axes E 2 and E 3 , which are depicted in Fig. 5 . To arrange the initial general geometry we rotate the system by means of the rotation tensor Q such that the unit basis vectors E i , i = 1 , 2 , 3 , become where Q = e 1 E 1 + e 2 E 2 + e 3 E 3 ,  (39) , L ± 0 are the directions of representative cross-links, which are rotationally symmetric with respect to ± E 1 . E R represents an arbitrary unit vector normal to E 1 . and hence with respect to spherical polar coordinates θ and φ shown in Fig. 5 e 1 = sin θ cos φE 1 + sin θ sin φE 2 + cos θ E 3 ,

ARTICLE IN PRESS
with e 1 now identified as the collagen fiber direction.
Here the directions of two symmetrically disposed families of cross-links are denoted by the unit vectors L + 0 and L − 0 , as distinct from the notation L ± used in (16) , so that where E R , which is an arbitrary vector orthogonal to E 1 , can be written as with φ 0 arbitrary. Then we define L ± according to L ± = Q L ± 0 = ± cos α 0 e 1 + sin α 0 e r , e r = Q E R = cos φ 0 e 2 + sin φ 0 e 3 .

(41)
Now on application of a deformation gradient F the invariant I 4 associated with the fiber direction is given by It follows from (41) 1 that F L ± = ± cos α 0 F e 1 + sin α 0 F e r .
Hence, the invariants I ± , and the quantities I ± 8 describing the coupling between the collagen fiber and cross-link directions are where for conciseness we have written s 0 = sin α 0 and c 0 = cos α 0 . Note that, in general, I + = I − and I + 8 = −I − 8 , which is unlike the case of the uniaxial tension considered in Section 3 .
Next we note the derivatives of the invariants I 4 , I ± and the quantities I ± 8 with respect to the right Cauchy-Green tensor C , i.e. JID: MPS [m3Gsc;August 16, 2019;16:12 ] ∂I ± ∂ C = c 2 0 e 1 e 1 ± s 0 c 0 ( e 1 e r + e r e 1 ) + s 2 0 e r e r ,

(48)
Now let us consider the strain-energy function (I 1 , I 4 , I + , I − , I + where we have used the abbreviations ψ 1 = ∂ /∂ I 1 , ψ 4 = ∂ /∂ I 4 , ψ I ± = ∂ /∂ I ± and ψ 8 ± = ∂ /∂ I ± 8 . This is the most general Cauchy stress expression for parallel collagen fibers with cross-links of the type indicated. For the related elasticity tensor in the material description see the Appendix. To recover the uniaxial case ( Section 3 ) from the general equations in this section we have e 1 = E 1 for uniaxial tension, and consequently F e 1 = λe 1 , F e r = λ −1 / 2 e r , F L ± = ±c 0 λe 1 + s 0 λ −1 / 2 e r .
which gives with (56) the same expression for σ as in (33) except that here the volume fractions are incorporated into the material constants μ, k 1 , ν and κ.

Planar formulation
Next we consider the situation in which the fibers and cross-links are restricted to the ( E 1 , E 2 ) plane and we define the fiber direction e 1 and its normal e r as e 1 = cos αE 1 + sin αE 2 , e r = − sin αE 1 + cos αE 2 , where α is the angle between the fiber direction and the E 1 axis (see Fig. 6 ).
With respect to e 1 and e r the cross-link directions L ± are defined by L ± = ±c 0 e 1 + s 0 e r , which has the same form as (41) 1 . The invariant I 4 = ( C e 1 ) · e 1 , as in (42) 2 , but with e 1 now defined by (58) . We also have F L ± = ±c 0 F e 1 + s 0 F e r , which is the same expression as (43) and the invariants I ± and the quantities I ± 8 are again given by (44) and (45) . The Cauchy stress tensor σ has the same form (49)  In particular, e 1 represents the fiber direction with unit normal e r with respect to background axes E 1 and E 2 , while e 1 makes an angle α with respect to the E 1 -direction. L ± represent the directions of two families of cross-links, and L ± make an angle α 0 with respect to the ± e 1 direction.

Simple shear case
For simple shear in the E 1 direction in the considered plane, the deformation gradient is given by F = I + γ E 1 E 2 , where γ is the amount of shear. It follows that F e 1 = e 1 + γ sin αE 1 , F e r = e r + γ cos αE 1 .

ARTICLE IN PRESS
Hence, from (68) we obtain In Fig. 7 we plot the dimensionless shear stress σ 12 = σ 12 /μ from (74) against the amount of shear γ in order to illustrate the dependence on the various parameters. Fig. 7 (a) shows how the results depend on the angle α 0 of the cross-links relative to the fiber direction for a fixed value of the fiber angle α and each of the other parameters, as specified in the caption. Fig. 7. Plots of the dimensionless Cauchy shear stress σ 12 = σ 12 /μ versus the amount of shear γ : (a) for three values of the cross-link angle α 0 ( π/12, π/6, π/4) compared with the plot for the case of no cross-links. On the basis of (74) the following parameters were used: k 1 = k 1 /μ = 1 , ν = ν/μ = 2 , κ = κ/μ = 1 , α = π/ 3 , k 2 = 0 . 1 ; (b) for four values of the dimensionless cross-link parameter ν (8.0, 5.0, 2.0, 0). On the basis of (74) the following parameters were used: k 1 = 1 , κ = 1 , α 0 = π/ 6 , α = π/ 3 , k 2 = 0 . 1 . JID: MPS [m3Gsc;August 16, 2019;16:12 ] Analogously to Fig. 4 , Fig. 7 (b) illustrates the dependence on the density of cross-links via the dimensionless parameter ν and fixed values of each of the other parameters. Once again it is clear that the response becomes stiffer as either the cross-link direction approaches the fiber direction or the density of cross-links increases.

Concluding remarks
In this paper we have developed two basic and preliminary continuum models that include the effects of collagen fiber cross-links on the elastic behavior of fibrous soft biological tissues. The first approach is a simple and purely phenomenological extension of a standard model that is used to describe the anisotropic behavior of fibrous tissues by including the dependence on the cross-link density in the anisotropic stiffness parameter of the standard model. This model has been examined for the case of homogeneous simple tension and it illustrates how the stiffness of the material increases with increases in the density of cross-links.
The second model is more general and aims to account for both the relative orientation and stiffnesses of (two families of) collagen fibers and cross-links, and their coupling using an invariant-based strain-energy function. The model is first applied to the case of simple tension in the fiber direction, which shows how the relative orientation of the cross-links and the fibers, and the cross-link stiffness affect the uniaxial response. Then, a more general 3D form with an arbitrary fiber direction and with two symmetrically disposed families of cross-links is provided, which is also suitable for finite element implementation, as needed for solving more complex 3D boundary-value problems. For the purpose of such an implementation, the elasticity tensor for the model in the material description is also provided.
The 2D (plane strain) specialization of the model is applied to the case of a simple shear deformation, and this shows again how the relative orientation of the cross-links and the fibers increases the overall stiffness as the relative orientation is reduced, as well as the cross-link stiffness.