New experimental results on collagen fibre dispersion in human arterial layers have shown that the dispersion in the tangential plane is more significant than that out of plane. A rotationally symmetric dispersion model is not able to capture this distinction. For this reason, we introduce a new non-symmetric dispersion model, based on the bivariate von Mises distribution, which is used to construct a new structure tensor. The latter is incorporated in a strain-energy function that accommodates both the mechanical and structural features of the material, extending our rotationally symmetric dispersion model (Gasser et al. 2006 J. R. Soc. Interface 3, 15–35. (doi:10.1098/rsif.2005.0073)). We provide specific ranges for the dispersion parameters and show how previous models can be deduced as special cases. We also provide explicit expressions for the stress and elasticity tensors in the Lagrangian description that are needed for a finite-element implementation. Material and structural parameters were obtained by fitting predictions of the model to experimental data obtained from human abdominal aortic adventitia. In a finite-element example, we analyse the influence of the fibre dispersion on the homogeneous biaxial mechanical response of aortic strips, and in a final example the non-homogeneous stress distribution is obtained for circumferential and axial strips under fixed extension. It has recently become apparent that this more general model is needed for describing the mechanical behaviour of a variety of fibrous tissues.
Physiological and pathological changes in the cardiovascular system directly influence the mechanical behaviour of arterial walls . It is, therefore, of crucial importance to improve understanding of the mechanical properties of the constituents of arterial walls, including the inherent features of anisotropy and nonlinearity. These properties, among others, pose formidable challenges in the constitutive modelling and numerical analysis of such tissues and can be clearly connected to the underlying structure of the tissues. The passive mechanical behaviour of an arterial wall is governed mainly by the matrix material (which consists of water, elastin, proteoglycans, etc.) and the collagen fibre reinforcement. The anisotropy is associated with the local mean alignment of the collagen fibres, which also stiffen their response when under tension, leading to their significant nonlinear characteristics. The fibres are not perfectly aligned but are dispersed around their mean direction, and the amount and character of the dispersion depends on the topography, the particular layer of the vessel considered and the respective (patho)physiological condition, inter alia. Fibre dispersion in arterial walls has been documented and analysed in, for example, [2–5]; for an overview of the structural quantification of collagen fibres in arterial walls, see . In particular, the study in  identified the presence of two fibre families in the intima, media and adventitia of human aortas; they are helically and almost symmetrically arranged with respect to the cylinder axis. Often a third and sometimes a fourth family is present in the intima in the respective axial and circumferential directions. Recent work  has revealed that, while helical fibre structures are present in human elastic arteries, in more muscular arteries (as for the murine basilar artery) and veins (such as the porcine jugular vein) a transition from the helical arrangement to two nearly orthogonal fibre families aligned in the circumferential and axial directions can be observed, and it is suggested that this is to ensure optimal efficiency of the vasculature. Observations of dispersion for other tissues, including the myocardium, corneas and articular cartilage, can be found in [8,9], [10,11] and , respectively.
Several mechanical models accounting for the dispersion of collagen fibres have been proposed. Fibre dispersion can be represented either directly by incorporation in a strain-energy function via a probability density function (PDF) or by a generalized structure tensor, for example. Following , these two approaches are referred to as ‘angular integration’ (AI) and ‘generalized structure tensor’ (GST). Cortes et al.  compared the results of the two formulations on the basis of the energy function introduced in  and the GST in . As was pointed out in, for example,  one approach is to consider the strain energy w(λ) of a single collagen fibre as a function of the fibre stretch λ and to integrate this over the unit sphere S to obtain the total free-energy function Ψf of the fibres per unit reference volume, i.e. 1.1where n is the number of fibres per unit reference volume, N is a unit vector describing the orientation of an individual fibre and ρ is the relative angular density of fibres normalized according to 1.2
The region of the unit sphere where fibres are in tension is defined by the set of N for which λ > 1, where λ = (C : N ⊗ N)1/2 is the stretch in the direction N, C is the right Cauchy–Green tensor and a double contraction is defined by C : N ⊗ N = (CN) · N. This is the AI approach.
In the GST approach, the energy function is associated with a GST H and is given by 1.3where H is defined by 1.4the latter following from (1.2).
By comparing these approaches, we see that (1.1) requires integration over the unit sphere at each point, while for the strain-energy function (1.3) no such integration is needed once H has been determined by (1.4)1 or prescribed. On the one hand, the AI approach allows identification and exclusion of individual fibres under compression, while, on the other hand, this is not quite so straightforward on the basis of (1.4)1. For a detailed discussion of the latter issue, see [17,18].
In chronological order, we now provide a short overview of the main existing continuum mechanical models which take fibre dispersion into account. Probably the study in  was the first to consider fibre dispersion in the analysis of fibrous connective tissues, with the tissue structure (fibre orientation) being accounted for in the strain energy via an orientation density function. The approach  incorporates a two-dimensional AI distribution in a strain-energy function based on the Beta distribution, but neglects the out-of-plane dispersion of the fibres, while a planar fibre dispersion was used in  with a Gaussian PDF to study the biaxial behaviour of arterial walls and aortic valves. In , the Gaussian distribution was used to compute the so-called splay invariants to represent two- and three-dimensional fibre dispersions, which the authors applied to aortic valve tissues.
In , the use of an exponential strain-energy function was motivated following , and this model was extended to the case of fibre dispersion in  based on the structure tensor (1.4)1. Therein, we used a rotationally symmetric distribution of collagen fibres. In addition, the model  was applied to several other tissues, including the cornea  and articular cartilage , while, in , the constitutive model of  was applied with the fibre orientations uniformly distributed over the azimuthal angle. A structure tensor based on a planar counterpart of that in  was defined in [27,28] and used in .
A PDF based on the von Mises distribution taking account of a non-rotationally symmetric fibre dispersion and based on the micro-sphere model was suggested in , and an ellipsoidal distribution with a power-law strain-energy function, based on the AI formulation, was employed in  and applied to cartilage. Based on the AI approach with the von Mises distribution Raghupathy & Barocas  derived a closed-form solution for a simple exponential fibre stress–strain law and applied their model to planar biaxial extension of a bioartificial tissue. The study of Federico & Gasser  adopted the model of Gasser et al.  and included the limiting case of an in-plane arrangement of fibres following the AI approach, which was also used in  with a planar von Mises distribution to examine the in-plane dispersion of collagen fibres. In , the Bingham distribution was used and was claimed to be ‘clearly superior to the π-periodic von Mises distribution in modelling the collagen organization in vascular tissue’, which we have not found to be the case in this study. The model  was extended in  to incorporate a higher order statistical measure of dispersion in order to reduce differences between the GST and AI formulations. In a recent paper , we have introduced a fibre dispersion model with a weighted energy function that allows the exclusion of fibres which are under compression. The model, based on the AI approach, has been developed for plane strain and for three-dimensional deformations appropriate for finite-element implementation (see also the discussion in  on the exclusion of compressed fibres using the GST approach). A summary of the main models discussed above is listed in table 1.
The model in  has proved to be very successful, but recently it was shown in  that an axisymmetric model of collagen fibre dispersion is not appropriate, and a more general dispersion model is required to accommodate the new findings. Hence, based on the structure tensor approach initiated in , we introduce here a bivariate von Mises distribution that enables the dispersion data to be captured. In particular, for arteries, the out-of-plane dispersion is relatively narrow while the in-plane dispersion is more significant. This work provides a natural extension of the constitutive setting documented in [14,15] to a more general context.
This work is structured as follows. In §2, we introduce the mathematical framework for describing fibre dispersion illustrated by the use of the bivariate von Mises distribution. This is particularly appropriate for use with the new experimental data, which indicate that the current models are not sufficiently general to capture the collagen fibre structure. Then we introduce a PDF which accounts for these experimental observations, and it is shown how the new model reproduces several special cases from the literature. The continuum mechanical framework associated with the dispersion structure introduced in §2 is presented in §3. In §4, the theory of §3 is applied in three representative numerical simulations with the aim of showing the efficacy and capability of the proposed structural model. Finally, §5 summarizes the proposed fibre dispersion model and discusses future developments of our work.
2. Mathematical representation of fibre dispersion
2.1. Structure tensor for the fibre dispersion
Motivated by the experimental results documented in , we introduce the coordinate system shown in figure 1, where the unit vector N is an arbitrary fibre direction in the reference (undeformed) configuration, expressed in terms of the two angles Φ and Θ by 2.1where , and e1, e2, e3 are unit rectangular Cartesian basis vectors. Locally, e1 and e2 define the tangential plane of a cylindrical coordinate system and e3 the corresponding outward radial direction. For a circular cylinder, e1 is taken to be the circumferential direction and e2 the axial direction. Although the coordinate system is similar to that introduced in , there is a subtle but important difference: our approach does not involve symmetry about a preferred direction since recent experimental results  have suggested that the fibre dispersion is not rotationally symmetric.
We describe ρ(N), which appears in (1.1), as the probability density of the fibre orientation N in the reference configuration as a function of Φ and Θ. In the following, we write either ρ(N) or ρ(Φ, Θ), and this is normalized according to (1.2), equivalently 2.2where dS = cosΘdΦdΘ. (Note that the usual spherical polar angles are π/2 − Θ and Φ.)
Another requirement is that the PDF has to be independent of the sense of N so that ρ(N) ≡ ρ(−N), which is equivalent to ρ(Φ, Θ) = ρ(Φ + π, −Θ). Based on the experimental results presented in , we introduce two additional symmetries, namely the in-plane symmetry ρ(Φ, Θ) = ρ(Φ + π, Θ) and the out-of-plane symmetry ρ(Φ, Θ) = ρ(Φ, −Θ).
It is assumed that the material behaviour does not depend on the sense of N, so that the strain energy depends on N only through the tensor product N ⊗ N (; see also [15,39]), via a symmetric second-order tensor, introduced in (1.4)1, which we now write as 2.3the components of which involve only ρ(N) and the sines and cosines of the angles Φ, Θ . The tensor H is a structure tensor involving the fibre dispersion via ρ(N). This is associated with a dispersion of fibres about a single mean direction, the mean direction of N, which, according to the data in , has a very small component out of plane. We therefore assume here that the mean fibre direction, say M, lies in the tangential plane defined by e1 and e2. An explicit expression for M will be written down below.
In terms of its components Hij, we write H = Hij ei ⊗ ej and note that, due to the symmetries of ρ(N), the off-diagonal components H13 and H23 vanish, and the only non-zero components of H are therefore H12 and the diagonal components H11, H22, H33. Thus, H has four non-zero components, and, because trH = 1 by (1.4)2, only three of them are independent: 2.4
In , it was also observed that the dispersions in the two planes are essentially independent, which means that the PDF can be decomposed in the form 2.5where ρip(Φ) and ρop(Θ) describe the in-plane and out-of-plane dispersions, respectively. Since the symmetry of the PDF has to be fulfilled, we require that
Without taking explicit account of these symmetries, the normalization in equation (2.2) is now written as 2.6
Without loss of generality, we can choose the normalization of ρip so that 2.7and hence (2.6) reduces to 2.8Let us now define 2.9which is a measure of the out-of-plane dispersion and is consistent with the definition, in slightly different notation, given in eqn (4.1) in  for the case in which the dispersion is rotationally symmetric. By definition, κop is non-negative and from (2.8) and (2.9) it must satisfy κop ≤ 1/2. Thus, 2.10If ρop ≡ 1, giving an isotropic distribution, then κop = 1/3. Note that if ρip = 1 then (2.7) is automatically satisfied, but on the other hand equation (2.7) does not necessarily imply that ρip = 1. However, when ρip = 1 we obtain H11 = H22 = κop, H33 = 1− 2κop and H12 = 0 and this corresponds to a transversely isotropic distribution with mean fibre direction e3. The lower limit κop = 0 corresponds to the case in which all the fibres are in the e3 direction (no dispersion), while the upper limit κop = 1/2 corresponds to a planar distribution in the (e1, e2) plane. For isotropy, we have ρip = ρop = 1 and κop = 1/3. In general, ρip and ρop are separate measures of the fibre dispersions, in plane and out of plane, respectively.
We now define and by 2.11and 2.12and hence by (2.11) 2.13Together, and characterize the in-plane dispersion. Note that when ρip = 1 and also when the mean in-plane direction is either e1 or e2, as discussed below.
In general, we now have 2.14
At this point, we note that the structure tensor H has the form 2.15and we emphasize, as noted above, that it is associated with a dispersion that has a single mean direction. We shall relate this to the model of dispersion discussed in  when considering special cases in §2.2.
An alternative way of writing the above representation is the form 2.16where the unit vector M = cosαe1 + sinαe2 is the in-plane mean fibre direction, α being the angle between M and e1, and Mn is a unit out-of-plane vector (see figure 2), while A and B are constants.
By comparing (2.15) and (2.16), we find that the components Hij are related to A, B and the angle α by 2.17from which we obtain 2.18and 2.19
By defining and , we obtain 2.20where and depend only on the in-plane distribution. Note that is invariant with respect to the rotation of the in-plane axes.
For any given distribution ρip, we can calculate and and hence the angle α and the constants and , while κop is obtained when ρop is given. Note that H contains only two independent parameters. Note also that for α = 0 or π/2 the dispersion parameter .
In §2.1.1, we consider a particular choice of the density functions ρip and ρop.
2.1.1. Using the bivariate von Mises distribution as the dispersion distribution
The work of Schriefl et al.  documented angular datasets for the in-plane collagen dispersions of the intima, media and adventitia of human non-atherosclerotic thoracic and abdominal aortas and common iliac arteries. The out-of-plane angles were measured separately for each layer, showing that the out-of-plane dispersions are very similar at all anatomic locations and for each layer. Moreover, each mean fibre angle was found to be very close to tangential. Motivated by these results we now model, as illustrative examples, each of ρip(Φ) and ρop(Θ) with π-periodic von Mises distributions , so the overall PDF is a bivariate von Mises distribution, featuring the symmetries that were discussed in §2.1.
For ρip(Φ), we consider the basic von Mises distribution 2.21where the so-called concentration parameter a is a constant and I0(a) is the modified Bessel function of the first kind of order 0 defined by 2.22which provides a normalization factor leading to (2.11) independently of a. Equation (2.22) corresponds to a distribution that is symmetric about Φ = 0. Plots of the circular distribution (2.25) for different concentration parameters a are shown in figure 3, and we note that, as a → ∞, ρip becomes a delta function. For this special case of symmetry, we use the notation κ11, κ22 and κ12. It follows from (2.25) using the definitions (2.15) and (2.16) that κ12 = 0 and 2.23where 2.24is the modified Bessel function of the first kind of order 1.
If the distribution is symmetrical about Φ = α instead of Φ = 0, then Φ is replaced by Φ − α in (2.21) and the appropriate values of and are given by 2.25 and is given by (2.13). From now on, whenever κ12 = 0 we use the notation κip instead of κ22.
Equation (2.8) is satisfied by taking ρop(Θ) to be a von Mises distribution of the form 2.26where b is a constant concentration parameter and erf is the error function defined by 2.27Note that (2.26) has a similar shape to (2.21) but is marginally different from the distribution used in , and gives a somewhat better fit to the narrow out-of-plane dispersion.
For the distribution (2.26), a closed-form expression for κop from (2.9) is obtained in the form 2.28
Figure 4 shows a plot of κop as a function of the concentration parameter b. Note in particular that κop = 1/3 when b = 0.
We use maximum-likelihood estimates to obtain the parameters a and b in the PDFs, ρip(Φ) and ρop(Θ), from the angular datasets of Schriefl et al. . Although it is possible to determine the parameters of a PDF by minimizing the sum of squared errors, this method has several disadvantages, as pointed out in , and hence we prefer to identify the parameters by using a maximum-likelihood estimate. In figure 5, we show experimental data from the adventitia of a human non-atherosclerotic abdominal aorta obtained from picrosirius polarization, in combination with a universal stage . To produce the fit, we used the maximum-likelihood estimate to obtain the concentration parameters a and b. In figure 5a, we show the in-plane bi-modal distribution of ρip of Φ with α = ±47.99° giving a = 2.54, which is obtained from (2.21) using Φ + α and Φ − α instead of Φ. In figure 5b, we show the out-of-plane probability density ρop of Θ obtained from (2.26), where b = 19.44.
2.2. Special cases of fibre dispersions
Our model includes several existing dispersion models as special cases of (2.16), which are discussed in the following section.
2.2.1. Transversely isotropic dispersion
In , we considered a transversely isotropic dispersion for which H has the form 2.29where κ is a single dispersion parameter associated with a transversely isotropic dispersion about M. Note that here we are taking M to lie in the (e1, e2) plane (figure 2). In , κ had the form where Θκ is measured from M and is different from the Θ used here. The connection between Θκ and Θ is obtained from M · N = cosΘκ, with N given by (2.1) and M = cosαe1 + sinαe2, in the form cosΘκ = cosΘ cos(Φ − α).
Equation (2.29) is recovered as a special case of (2.16) by taking κ = 1 − 2κop, which corresponds to A = κ, B = 1 − 3κ.
2.2.2. Perfect alignment
If both concentration parameters a and b become infinite, there is no dispersion in either plane and we obtain the model proposed in . With a → ∞, ρip becomes an in-plane delta function, and with b → ∞ we have κop → 1/2. The structure tensor is then H = M ⊗ M, and all fibres are oriented in the in-plane direction of M. This corresponds to A = 0, B = 1 in (2.16).
2.2.3. Isotropic dispersion
An isotropic fibre dispersion is represented by a uniform dispersion in each plane, meaning that ρ(Φ, Θ) is independent of Φ and Θ, with ρip = ρop = 1, so that, for the von Mises distributions, a = b = 0, while κop = 1/3, and the structure tensor is given by H = (1/3)I. Thus, there is no preferred direction (κ = 1/3), with A = 1/3, B = 0 in (2.16).
2.2.4. Planar dispersion
A dispersion with fibres oriented only in plane was presented in [27,39]. In this case, there is no out-of-plane contribution to H, which can be written as 2.30where 1 is the two-dimensional identity in the considered plane with the normalization and dispersion parameter given by 2.31Note that the normalization (2.31)1 is equivalent to (2.7) bearing in mind the different ranges of angles used, and κ is equivalent to as defined in (2.11)2. The in-plane PDF ρ(Θ) satisfies ρ(−Θ) = ρ(Θ). For the von Mises distribution (2.26), this corresponds to b → ∞ and κop → 1/2. Equation (2.30) is obtained from (2.16) by setting , B = 1−2κ. Note that this is a two-dimensional distribution but its application is not restricted to use in two dimensions.
2.2.5. Planar isotropic dispersion
If a dispersion features perfect out-of-plane alignment, b → ∞ and κop → 1/2, and is fully dispersed in plane so that a → 0, ρip = 1, then it is planar isotropic. For this case, the structure tensor is simply H = (1/2)1, corresponding to A = 1/2, B = 0 in (2.16).
2.2.6. Summary of the special cases
The special cases discussed above are summarized in table 2. Figure 6 is a visualization of ρ(N)N for (a) the general case for which H is given by equation (2.16), (b) the transversely isotropic dispersion given in §2.2.1 with H given by (2.29), (c) the case of perfect alignment according to §2.2.2, (d) the isotropic case according to §2.2.3 and (e) the case of planar dispersion given by equation (2.30).
3. Continuum mechanical framework
We consider a (stress-free) reference configuration, denoted Ω0, and a deformed (spatial) configuration, denoted Ω. The deformation map χ(X) transforms a material point X ∈ Ω0 into a spatial point x∈Ω. With this deformation map we define the deformation gradient F(X) = ∂χ(x)/∂X and its determinant J = detF(X), where J is the local volume ratio; we require J > 0.
Following [42,43], we apply the multiplicative decomposition of F into a spherical (dilatational) part J1/3I and a unimodular (distortional) part , with . We define the right Cauchy–Green tensor and its modified counterpart as C = FTF and , respectively, with the related invariants I1 = trC and .
Since we treat an artery as an elastic material, we assume the existence of a strain-energy function Ψ(C, H4, H6) that depends on the macroscopic deformation through C and the underlying tissue structure through the structure tensors H4 and H6, which describe the fibre alignment and dispersion for two fibre families. Based on (2.16), these are defined by 3.1and M4 and M6 lie in the (e1, e2) plane, and Mn normal to that plane.
For computational purposes, we assume that it is possible to split the strain-energy function into two parts as , as shown in [5,44]. The function Ψvol is a purely volumetric contribution while Ψiso represents the energy contribution of an isochoric (volume preserving) deformation through . The second Piola–Kirchhoff stress tensor S is given by S = 2∂Ψ/∂C. Using the decoupled form of Ψ, we can identify two stress contributions so that S = Svol + Siso. Using well-known results from tensor analysis (e.g. ), and the chain rule, we obtain 3.2where p = dΨvol(J)/dJ is the constitutive equation for the hydrostatic pressure and 3.3is the so-called ‘fictitious’ isochoric second Piola–Kirchhoff stress tensor. The deviator in the Lagrangian configuration is defined by , where is a projection tensor that furnishes the correct deviatoric operator in the Lagrangian setting, and is a fourth-order identity tensor.
The related elasticity tensor in the Lagrangian description is written in the decoupled form 3.4
According to , we have the specification 3.5and 3.6In (3.5), we have used the notation , where the symbol denotes the tensor product according to the rule 3.7and the scalar function is defined by with the constitutive equation for p given in the line after equation (3.2). In (3.6), we have also used the definitions 3.8where is the fourth-order fictitious elasticity tensor, is the trace and is the modified projection tensor of fourth order. The related spatial stress and elasticity tensors may be derived by push-forward operations on (3.2) and (3.4)1, respectively . It should be noted here that the compressible formulation is introduced for computational purposes, and the incompressibility condition has to be enforced by a numerical scheme, one example of which is the augmented Lagrangian method (e.g. ).
We emphasize that for a specific material we need to specify Ψiso and hence to calculate its derivatives with respect to , which affects the two expressions (3.3) and (3.8)1. In §3.1, such a specification is provided.
3.1. Anisotropic strain-energy function
Each of the structure tensors H4 and H6 depends on two dispersion parameters. We assume that these are the same for each fibre family. Our approach follows the work of Holzapfel et al.  in which the contributions Ψg and Ψf of the ground matrix (non-collagenous material) and the fibres to the strain energy are added. The artery is treated as an incompressible, elastic and fibre-reinforced material with the fibre dispersion accounted for both in plane and out of plane. Hence, superposition of energies reads 3.9
Following [14,46], we model the ground substance with a neo-Hookean material , where the parameter c is the shear modulus in the reference configuration. For the fibre contributions Ψfi, we adopt the exponential functions  3.10where k1 is a parameter with the dimensions of stress and k2 a dimensionless parameter, while 3.11is a Green–Lagrange strain-like quantity which can be interpreted as an averaged or weighted fibre strain, depending on the fibre dispersion through the structure tensor Hi and the (isochoric) macroscopic deformation through .
Since trHi = 1, we can write . Using the definitions (3.1) of the structure tensors H4 and H6, we obtain 3.12where 3.13The invariants and are the squares of the stretches in the directions Mi and Mn, respectively. A summary of the parameters used is provided in table 3.
The strain-energy function used in the proposed model reads 3.14
Following , we make the common assumption that the fibres do not resist any compression and are only active in tension. The invariants Ii are used as switches between fibre compression and tension so that the anisotropic part Ψfi only contributes to the strain energy when I4 > 1 or I6 > 1. If one or more of these conditions is not satisfied then the relevant part of the anisotropic function is omitted from (3.14). For example, if I4 and I6 are less than (or equal to) 1, then the tissue response is purely isotropic. For discussion of subtle points regarding the choice of switching criteria, see .
In the expressions for the stress and the elasticity tensor, we need to calculate and . By using the chain rule, these are obtained as 3.15where Hi is given by (3.1), and we have introduced the notations 3.16
4. Representative examples
Having identified the dispersion parameters of the model in §3, we now identify the mechanical parameters by fitting the constitutive model to uniaxial data. The proposed constitutive model is then implemented in the finite-element analysis program FEAP, and the parameters obtained are used to obtain results for a simple biaxial extension of an eight-element cube and compared with the analytical (Matlab) result. Finally, numerical results are obtained for the stress distribution in the non-homogeneous extension of strips of an adventitial layer cut out along the axial and circumferential directions.
4.1. Parameter fitting to experimental data
In this example, we consider the purely incompressible formulation where the strain-energy function is characterized by Ψ = Ψ(C, H4, H6). The second Piola–Kirchhoff stress tensor S is then given by 4.1where p in this case denotes the Lagrange multiplier required to enforce incompressibility. From (4.1), the Cauchy stress tensor can be computed simply by σ = FSFT.
The strain-energy function in the present formulation reads 4.2where 4.3The anisotropic term in (4.2) only contributes when the fibres are extended, i.e. when I4 > 1 or I6 > 1. For the case that one or more of these conditions is not satisfied then the relevant part is omitted from (4.2).
We now consider a tissue specimen with two fibre families in the reference configuration with mean directions illustrated in figure 7. The unit vectors M4 and M6, which are symmetrically disposed in the circumferential/axial plane, each make an angle α with the circumferential direction, so that 4.4while the normal direction to the plane is [Mn] = [0, 0, 1]T. The deformation gradient matrix and the right Cauchy–Green matrix are 4.5where λ1, λ2, λ3 denote the principal stretches. Hence, the required invariants read 4.6 4.7 4.8Since Ii is the same for i = 4 and 6, we obtain , where here, in analogy with (3.16)1, , i = 4, 6. By using Ei in the strain-energy function (4.2), we obtain the derivatives of the strain-energy function as 4.9Hence, the non-zero components of the Cauchy stress are 4.10 4.11 4.12The constants A and B can be deduced from (2.17) and (2.18) as 4.13
Together with the incompressibility condition (λ1λ2λ3 = 1), the implicit equations (4.10)–(4.12) with σ22 = σ33 = 0 can then be used to obtain p, λ2 and λ3 in terms of λ1, thus giving an expression for σ11 in terms of λ1, the material parameters c, k1 and k2 and the structural parameters κip, κop and α.
To determine the material parameters of our model, we use the structural data provided in figure 5 (mean fibre angle α and concentration parameters a, b) and experimental data (from uniaxial tension tests) from the adventitia of a human non-atherosclerotic abdominal aorta. A least-squares objective function is chosen as the L2-norm of the error between the model prediction of the Cauchy stress in the circumferential and the axial directions, and the corresponding experimental data. We use the function lsqnonlin in Matlab (v. R2010b; The MathWorks Inc., Natick, MA, USA), which uses the Levenberg–Marquardt algorithm to find the minimum. The dispersion parameters κip from (2.23)2 with the value of a, and κop from (2.28) with the value of b, can be determined from the imaging data provided in figure 5. Imaging and material data were obtained from different but comparable tissues.
To quantify the goodness of the fit, we calculate the coefficient of determination R2, which is given by 4.14where Serr and Stot are the sums of squares of the differences between model/experiment and the mean of experiment/experiment, respectively . We obtained a good fit with R2 = 0.998. The parameters identified are listed in table 4. Figure 8 shows the result of fitting the proposed model to data from uniaxial tension tests.
4.2. Biaxial extension: comparison of Matlab and finite-element results
In this section, we demonstrate the efficacy of the proposed constitutive model by using eight hexahedral Q1/P0 finite elements for a cube (10 × 10 × 10), which is reinforced by two symmetric fibre families, as depicted in figure 7. The cube is subjected to homogeneous biaxial extension in the 1,2-plane, which here we call ‘in plane’, as shown in figure 9; the 2,3-plane we call ‘out of plane’. The mesh is unstructured according to figure 9. The material parameters used are those documented in table 4. We investigate four different cases of fibre dispersion: (i) structural parameters taken from table 4, which is the reference case; (ii) high alignment out of plane and isotropy in plane; (iii) less alignment out of plane and isotropy in plane; (iv) isotropy in both out of plane and in plane. For a summary see table 5.
The maximum of each stretch, (λ1, λ2) in the 1 and 2 directions, respectively, is 1.15 and the corresponding Cauchy stresses (σ11, σ22) are computed from equations (4.10)–(4.12) using Matlab. In addition, these relationships are compared with results obtained from a finite-element computation using FEAP . The cube is stretched simultaneously in the 1 and 2 directions under displacement-driven conditions using the Newton–Raphson method. Figure 10a,b shows plots of the Cauchy stresses σ11 and σ22 versus the stretches λ1 and λ2, respectively. Note that, for case (iv), the stresses in the 1 and 2 directions are the same since the in-plane and out-of-plane dispersions are isotropic.
4.3. Extension of adventitial strips
In this section, we illustrate the results of the finite-element implementation of the proposed constitutive model, simulating uniaxial extension tests related to experiments on strips performed in our laboratory in Graz, Austria. The strips were taken from an adventitial layer cut out along the axial and circumferential direction of a human non-atherosclerotic abdominal aorta. The strips, each with initial length, width and thickness of 10.0, 3.0 and 0.5 mm, respectively, were subjected to a stretch of 1.3. Both ends of each strip were constrained so as to model the mounting in the testing machine and were not allowed to deform. The resulting deformation of each strip was therefore non-homogeneous. We assume uniform material parameters over the adventitial strips with values provided in table 4. Two symmetric fibre families, as shown in figure 7, are assumed to make an angle α of ±47.99° with the circumferential direction and to show a distribution characterized by κip = 0.116 and κop = 0.493, as also provided in table 4. We use 3200 hexahedral elements, applying the mixed Q1/P0 element throughout the simulation. Figure 11 shows the element results of the circumferential and axial specimens, subjected to a stretch of 1.3. Both circumferential and axial specimens show that the Cauchy stresses in the adventitia can be modelled within the experimentally predicted range of stresses (compare with the experimental results in figure 8). It turns out that the (axial and circumferential) stresses in the middle of the specimens are slightly higher than the model stresses provided in figure 8, which is due to the boundary conditions used in the present example.
As several previous approaches have shown (e.g. [15,25,37]) incorporation of fibre dispersion into a continuum mechanical framework for soft biological tissues is a challenging and important task. Indeed the mechanical response of such tissues depends significantly on the tissue structure, in particular the arrangement of the fibre dispersion.
One of the main goals of our previous papers, e.g. [14,15], has been to incorporate the structure of biological tissues into our models in order to capture the physiological and pathological mechanical mechanisms. Recent experimental data have shown the need for a more general model that takes account of the non-symmetric arrangement of collagen fibres. This motivates the introduction, in this paper, of the bivariate von Mises distribution for describing the collagen fibre dispersion. By fitting the bivariate PDF to angular distribution data gained from imaging/histological analysis, we can determine structural parameters that can be integrated within the framework of continuum mechanics, in particular hyperelasticity. This allows us to define a strain-energy function from which the stress and the elasticity tensors can be computed, thus facilitating an efficient implementation of the model into a finite-element code.
The proposed constitutive model introduces a new structure tensor which incorporates in-plane and out-of-plane fibre dispersions in a clear and simple way, leading to an invariant-based formulation of the strain-energy function, generalizing our previous models [14,15]. The strain-energy function can incorporate different structure tensors for different fibre families. In particular, the work of Schriefl et al.  showed that the number of fibre families depends on the location of the artery and the type of layer, although in most cases two collagen fibre families were reported. For the most general case considered here the structure tensor has three independent components, and these reduce to two when there is in-plane symmetry (with parameters κip and κop) or one when there is rotational symmetry about a single preferred direction (with parameter κ). Various special cases have been highlighted in §2.2, capturing the cases of rotational symmetry, in-plane dispersion, transverse isotropy and isotropy.
We have constructed a specific energy function that illustrates the efficacy of the model using three specific examples. In the first example, we have demonstrated a good fit of the model to experimental data using identified structural parameters as a fixed set. The model was then used to illustrate the effect of fibre dispersion on biaxial extension experiments, including a comparison of analytical and finite-element results. Finally, we have analysed the three-dimensional non-homogeneous stress response for a given overall extension of adventitial strips in the circumferential and axial directions obtained from a human abdominal aorta. The resulting stress responses deviate significantly from each other due to the different mechanical and structural properties in the two directions.
More data are required to determine the detailed dispersion of collagen fibres, not only for the arterial wall but also for other types of fibrous tissues. To inform the modelling process, second-harmonic generation, for example, in combination with optical clearing  is a powerful technique for obtaining collagen fibre dispersion data from various types of tissues.
In order to analyse the data within a continuum framework, two main approaches are used. First the AI approach, which is computationally expensive because it involves an integration over the unit sphere. It does, however, allow the exclusion of fibres undergoing compression. With the GST approach, it is possible to extract structural parameters and to include them in the continuum mechanical framework to account for the fibre dispersion. For an artery that exhibits a dispersion with a strong alignment in plane, the GST seems to be a feasible and very efficient method, although it does require a minor modification in order to exclude fibres which are under compression .
For further development, as mentioned above, many more data are needed to inform the modelling process. In particular, there is a pressing need to obtain in vivo data in order to construct more realistic models of tissue and organ mechanics. The model can also be extended to incorporate inelastic effects such as damage, viscoelasticity and muscle activation.
- Received March 3, 2015.
- Accepted March 23, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.