Structurally motivated material models may provide increased insights into the underlying mechanics and physics of arteries under physiological loading conditions. We propose a multiscale model for arterial tissue capturing three different scales (i) a single collagen fibre; (ii) bundle of collagen fibres; and (iii) collagen network within the tissue. The waviness of collagen fibres is introduced by a probability density function for the recruitment stretch at which the fibre starts to bear load. The three-dimensional distribution of the collagen fibres is described by an orientation distribution function using the bivariate von Mises distribution, and fitted to experimental data. The strain energy for the tissue is decomposed additively into a part related to the matrix material and a part for the collagen fibres. Volume fractions account for the matrix/fibre constituents. The proposed model only uses two parameters namely a shear modulus of the matrix material and a (stiffness) parameter related to a single collagen fibre. A fit of the multiscale model to representative experimental data obtained from the individual layers of a human thoracic aorta shows that the proposed model is able to adequately capture the nonlinear and anisotropic behaviour of the aortic layers.
Integrating structural information to mechanical constitutive models of arterial tissues can help us to better understand its mechanical behaviour in physiological and pathological conditions. The first models published to describe the mechanical behaviour of arteries are purely phenomenological . The widely used model by Holzapfel et al.  includes structural information by considering the contribution of the collagen fibres in two distinct directions in space. Gasser et al.  add a von Mises distribution to describe a rotationally symmetric dispersion of the collagen fibres, whereas Alastrué et al.  use a micro-sphere approach to model blood vessels comparing a phenomenological strain–energy function and a worm-like chain model for the collagen fibres. Federico & Gasser  also account for the fibre distribution by defining a probability density function (PDF) on the unit sphere for a phenomenological strain–energy function of the collagen fibres. Gasser et al.  describe the fibre orientation on the unit sphere with a Bingham distribution to distinguish the in plane from the out of plane fibre distribution, whereas Holzapfel et al.  use the bivariate von Mises distribution to capture non-symmetric dispersion and to construct a new structure tensor. Zulliger et al.  explicitly define fractions of elastin and collagen in their model. They account for the waviness of collagen fibres folding the strain–energy function of a single fibre with the fibre distribution. The model by Cacho et al.  also accounts for the waviness of the collagen fibres by defining a distribution function for the recruitment stretch. However, both models [8,9] do not account for the three-dimensional fibre distribution in the tissue.
Collagen triple helices, i.e. collagen molecules, consist of three monomers and are cross-linked and assembled into fibrils. These fibrils constitute fibres. Tensile experiments with collagen I are present in the literature on all levels, from monomers to fibres [10–13]. Miyazaki & Hayashi  tested fibres which they extracted from the patellar tendon of a Japanese white rabbit. These fibres are stretched with glass microtubes. The nominal stress–strain curves soften with increasing strain suggesting a linear relationship between the Cauchy stress and the stretch. Sopakayang et al.  present data for fibres reconstituted with collagen and obtained from the tail tendon of Sprague–Dawley rats. Their results are similar to those in .
In this study, we propose a constitutive model that includes the waviness of collagen fibres and accounts for their three-dimensional distribution. Hence, two functions are introduced to describe the probability density ρr for the recruitment stretch, and the orientation distribution ρ of two families of fibres provided by the bivariate von Mises distribution. The constitutive model only uses two material parameters, namely a shear modulus μm of the matrix material and a (stiffness) parameter μf related to a single collagen fibre. Structural information is included on three different scales (i) the smallest scale is a single collagen fibre that only bears load when the stretch is larger than the recruitment stretch; (ii) the intermediate scale describes the bundle of collagen fibres by folding the strain energy of the collagen fibre with the distribution function for the recruitment stretch; and (iii) the largest scale describes the mechanical behaviour of the collagen at the tissue level where the orientation distribution of the collagen fibres is defined. An additive decomposition of the strain energies for the matrix material and the collagen fibres yields the strain–energy function for the tissue. The fraction of the constituents are accounted for by volume fractions for the matrix material and the collagen fibre bundles. Compared with the widely used models [2,3], the proposed model depends on more experimental data and the effort for numerical implementation is higher. However, the microstructure and the fibre distribution in particular are described more accurately.
In addition, we investigate the mechanisms on the smaller scales resulting in the characteristic tissue-level mechanical response. This involves a discussion of the experimental results documented by Schriefl et al.  suggesting that the media of the human iliac artery contains only one family of fibres. We show that the measured distribution can be modelled with a fibre distribution function containing two families of fibres with distinct mean fibre directions.
We hypothesize that by building a material model with the described hierarchical structural information we may predict the macroscopic behaviour of some soft biological tissues, in general, and of arterial tissues, in particular, when model parameters are identified from experiments. We confirm our hypothesis by comparing the proposed model with experimental data from uniaxial tensile tests of human arterial layers, i.e. intima, media and adventitia.
2. Material and methods
We investigate a material under deformation with fibres embedded in an isotropic matrix material and the fibres deform affinely with the matrix. We describe the transformation from the reference configuration to the current configuration with the deformation gradient F. For an incompressible deformation, the volume ratio J = det F = 1 holds. Subsequently, we use the right and the left Cauchy–Green tensors, i.e. 2.1
The first invariant of these tensors is defined as their trace, i.e. I1 = trC = trb. Consider now an arbitrary direction vector M in the reference configuration representing the orientation of a fibre. The mapping to the current configuration is then m = FM. The square of the stretch λ in the considered direction is 2.2Hence, a fibre initially oriented in the direction M is elongated by a factor of λ and oriented in the direction m.
2.2. Strain–energy function
The given model has three different scales, as depicted in figure 1, namely the single (wavy) collagen fibre, the bundle of collagen fibres and the collagen network within the tissue, i.e. the tissue with two fibre families and a surrounding non-collagenous matrix material. The related strain energies are described in the following.
2.2.1. Strain energy of a single (wavy) collagen fibre
The collagen fibre is on the smallest scale, see figure 1a, and is crimped in the reference configuration. It does not bear any load in compression and therefore only contributes to the strain–energy function when stretched after it is uncrimped. When the fibre is uncrimped, we say it is recruited, and the stretch of the tissue at that time is called the recruitment stretch λr of this particular fibre. The force–stretch relationship of a single fibre is assumed to be linear [14,15], and thus the related strain energy wf of the fibre follows as 2.3where μf is a material parameter and λ is the stretch of the tissue. The fibres are in compression for λ < λr and then wf = 0.
2.2.2. Strain energy of a collagen fibre bundle
A collagen bundle is composed of several individual fibres with a varying degree of crimp, see figure 1b, and hence a distribution of recruitment stretches. The strain energy wb of the fibre bundle is 2.4where ρr is the PDF for the recruitment stretch λr. This relation is equivalent to the convolution wb(λ) = wf(λ, 1) * ρr(λ).
22.214.171.124. Recruitment stretch distribution
We choose the beta distribution 2.5to describe the PDF for the recruitment stretch λr, where β(a1, a2) is the beta function and a1 and a2 are parameters to control the shape of the beta distribution. Equation (2.5) is defined on the interval 0 < λr ≤ 1. We transform it to the interval λ1 ≤ λr ≤ λ2, where λ1 is the stretch at which the first fibre is recruited and λ2 is the stretch at which the last fibre is recruited. Thus, 2.6Advantages are that the beta function is defined on a closed interval, that the shape can easily be controlled using the two parameters a1 and a2 and that symmetric and non-symmetric distributions can be described.
2.2.3. Strain energy of the tissue
A collagen fibre in the tissue is oriented in space and characterized by the vector M. The relative density of the fibres in a specific direction is described by the orientation distribution function (ODF) ρ(M). We require that the ODF fulfils the normalization condition 2.7where Ω is the surface of a unit sphere. Note that with the given normalization we have ρ = 1 for isotropy.
The total strain energy per unit reference volume for all collagen fibres in the tissue is  2.8where n is the number of fibre bundles per unit reference volume and is, therefore, closely related to the volume fraction ν of the fibre bundles. Consider the volume of a collagen fibre bundle to be V0, then ν = nV0 holds. The strain–energy function for the tissue composite is an additive decomposition of the strain–energy functions of the matrix material Ψm and the collagen fibres Ψf, i.e. 2.9Using the neo-Hookean constitutive model for the matrix, one has  2.10where (1 – ν) and μm are the volume fraction and the shear modulus of the matrix material, respectively. Because I1 and λ are both functions of the right Cauchy–Green tensor C, the strain–energy function (2.10) depends on the deformation through C.
126.96.36.199. Orientation distribution function of fibres
The orientation M of a collagen fibre can be described in a spherical coordinate system with the polar angle φ∈[0, 2π) and the azimuth angle , corresponding to the in plane and out of plane angles, respectively. Measurements of the fibre orientation [16,18] suggest two fibre families in arterial tissue with a non-symmetric distribution around a mean fibre direction. This mean fibre direction is identified by α which is the angle between the circumferential direction of the arterial wall and the direction with the highest density of fibres, the ‘mean fibre angle’ (figure 1c). A suitable function to describe the orientation distribution of one fibre family, i.e. ρ1, is the bivariate von Mises distribution according to 2.11where c1 is a normalization factor. The von Mises distribution is aligned such that b1 describes the in plane distribution of the fibres and b2 the out of plane distribution of the fibres.
For two families of fibres, we superpose the bivariate von Mises distributions ρ1 with one mirrored around the circumferential direction ρ2 to obtain ρ = ρ1 + ρ2, i.e. where the normalization factor c = 2c1 accounts for the additive decomposition of the two distributions such that the integral (2.7) holds.
2.3. Cauchy stress tensor
We consider the arterial wall as an incompressible material. Hence, the Cauchy stress tensor σ is obtained by 2.13where p is a Lagrange multiplier associated with the incompressibility constraint which can be interpreted as a hydrostatic pressure, and I is the second-order identity tensor. The extra stress tensor is calculated as, e.g. , 2.14where the prime denotes the derivative with respect to λ. While the first term in (2.14) describes the behaviour of the matrix material, the second term is associated with the contribution of the collagen fibres.
2.3.1. Numerical treatment of the integrals
To solve the integrals in (2.7), (2.14) over the unit sphere, we use the efficient method suggested by Bažant & Oh , with m = 42 distinct direction vectors Mi, i = 1, 2, … ,m, and with different integration weights qi. The integral over some tensor-valued function A(M) can then be approximated with a sum 2.15
Owing to the symmetry of the integration scheme, we need only 21 weights when multiplying the weights qi by two. The direction vectors and corresponding weights are given in .
2.3.2. Fit to experimental data
We use experimental data of the three layers of large (human) arteries to determine the structural parameters of the model. First, we fit the bivariate von Mises function to the data of the collagen dispersion, and then we identify the parameters of the Cauchy stress tensor with the data from a uniaxial extension test. Both fits are performed by means of the ‘trust-region reflective’ algorithm in the MATLAB's1 lsqnonlin function. The goodness-of-fit is evaluated with the parameter of determination R2.
We mimic the deformation mode of the experiments (uniaxial extension) and denote the stretch in the direction of the applied force λu. Considering incompressibility, the deformation gradient is then given as . The transversal stretch λt needs to be numerically determined through (2.13)1 and (2.14) considering zero stress boundary conditions normal to the applied force direction.
2.4. Model sensitivity
After fitting to experimental data, we want to examine the sensitivity of the proposed model to changes in its parameters. We adopt the widely used local method, see, e.g. , to estimate the influence of small changes from the identified parameter values on the Cauchy stress. To this end, we calculate the partial derivatives of the Cauchy stress with respect to the nine parameters pj (j = 1, 2, … , 9) and normalize by the absolute values to obtain dimensionless sensitivities 2.16for every stretch λu, where i = circ, ax. The partial derivatives are estimated with the method of finite differences.
3.1. Parameter study for the recruitment stretch distribution ρr
Here, we study the influence of the PDF ρr for the recruitment stretch, i.e. equation (2.6), on the mechanical behaviour of the tissue. In order to illustrate the effects, we choose an equi-biaxial tension mode of an incompressible material such that , where λe is the (equi-biaxial) stretch in the two directions. The parameters a1 and a2 in (2.6) define the symmetry/non-symmetry of the function and the interval λi = λ2 − λ1 varies the ‘slenderness’ of the function. All examples start with the two material parameters μm = 0.01 MPa, μf = 1 mN μm and the structural parameters n = 10−3 μm−3, λ1 = 1, λi = 0.5, α = π/6, b1 = 1 and b2 = 5. The Cauchy stress components in the circumferential σcirc and axial σax directions versus the (equi-biaxial) stretch λe are reported.
By increasing the parameters a1, a2 simultaneously, whereby we have chosen a1 = a2, the ‘slenderness’ of the PDF ρr for the recruitment stretch λr is significantly increased (figure 2a). Figure 2b shows the associated Cauchy stress evolution. Note that with a1 = a2 = 1 one obtains an equal distribution. The increase in the slenderness has a notable effect in the low stress region; a slender beta distribution results in a low initial stiffness. The stiffness increases faster resulting in a more J-shaped curve such that the stiffness and the stress values are the same for stretches when all fibres are recruited, i.e. for λe > 1.5.
Figure 2c,d shows the effect that occurs by increasing the interval λi from 0.2 to 0.5, and by taking a1 = a2 = 2. The shapes of the stress responses result in equal slopes after the last fibre is recruited; however, the stress values for a given stretch are higher for smaller intervals λi.
By increasing either a1 or a2 from 2 to 8 and by keeping the other parameter constant at 2, we obtain a non-symmetry in the PDF ρr, see figure 2e, whereas figure 2f shows the effect on the stress–stretch behaviour for varying a2. The impact is similar to a combination of changing the slenderness and changing the interval. The stiffness in the low stress region is affected and the stress values for a given stretch change. As in the above-discussed cases, the stiffness is the same after all fibres are recruited.
3.2. Parameter study for the orientation distribution function ρ of fibres
Now, we study the influence of the fibre ODF ρ, i.e. equation (2.12), on the stress–stretch behaviour of the tissue. Specifically, we investigate the parameter for the in plane distribution of the collagen fibres which is most relevant in a physiological sense. The material and structural parameters are the same as described in §3.1. The in plane distribution is determined by the parameter b1 ≥ 0, where b1 = 0 describes an isotropic in plane distribution, and is a slender distribution around the principal fibre direction with the mean fibre angle α. By fixing b2 at 5, three different values of b1 are compared in figure 3. Figure 3a,b shows profiles of the in plane and the out of plane distributions, respectively. The Cauchy stress for the equi-biaxial extension test is shown in figure 3c. The stress–stretch response changes strongly from an isotropic response to a highly anisotropic one, for an increase in b1.
3.3. Fit to experimental data
Here, we aim to show the reliability and capability of our proposed multiscale model. Therefore, we use experimental data from human thoracic aortas, as documented by Weisbecker et al. , and describe the mechanical response. The data in reference  were collected to calibrate a constitutive model for the softening behaviour of arterial tissues. For the present analysis, we just use the data of the primary loading curve; as a representative example we use the data from donor IX.
Figure 4 shows the comparison between the experimental data and the multiscale model. The model is able to adequately capture the nonlinear and anisotropic behaviour of the aortic tissue for all three layers. Table 1 summarizes all model parameters and indicates which parameters are obtained from the literature and which from the model fit. We group the parameters into three categories. The material parameters are followed by the structural parameters, i.e. the recruitment stretch distribution and the orientation distribution of fibres.
The initial shear modulus of the isolated (non-collagenous) matrix material from the thoracic aortic media is reported in Weisbecker et al. . In these samples, the collagen was digested by enzymatic treatment, hence this value represents the shear modulus of the volume fraction of the non-collagenous tissue. Therefore, in our model, the parameter corresponds to (1 – ν)μm. To the authors' best knowledge, data for the intima and the adventitia for the human aorta are not yet available in the literature. Thus, we use the same value for all three layers and that corresponds to the median value of the seven samples documented in .
We determine the value for μf from the experimental data documented in . The authors report a mean Young's modulus of E = 54.3 MPa for the individual collagen fibres. The relationship between E and the spring constant k is EA/L, where A is the fibre cross section and L is the fibre length. Now, we describe the elongation of the fibre by means of the stretch λ instead of l − L, where l is the length of the fibre in the current configuration. In order to obtain an equivalent expression for the strain energy, the spring constant is converted to μf = kL2 and thus μf = EAL. According to Miyazaki & Hayashi , the mean diameter of the fibres is D = 1.01 µm, and the mean reference length is L = 120 µm. With the modulus reported above, we obtain μf = 5.22 × 10−6 N mm = 5.22 mN µm, where the reference cross section is A = D2π/4. The intima and adventitia predominantly contains collagen type I, whereas the medium contains predominantly collagen type III . The literature reports data only for collagen type I, therefore, here, we use these data for all three layers.
The fitting parameter n describes the number of fibre bundles per reference volume of the tissue. The values in table 1 suggest that the diameter of the fibre bundles is considerably larger than 1 µm. This is a plausible result when comparing with, e.g. figure 5 in , where collagen fibre bundles of the intima are shown.
The lower bound of the recruitment stretch is described by λ1. We set λ1 = 1 by assuming that some fibres may be recruited immediately after initiation of loading. The interval of the recruitment stretch is defined by the parameter λi = λ2 − λ1. An automated method for the measurement of the fibre crimp distribution by analysing histological cuts was proposed earlier . Elbischger et al.  obtained a recruitment stretch range of 1 ≤ λr ≤ 1.5 for the adventitia of human iliac arteries.
Regarding the parameters of the beta distribution, we assume that a1 = a2. Comparing the stress–stretch response in figure 2, we see that a change in the interval λi has a similar effect on the mechanical behaviour as a change of the ratio a2/a1, namely an increase in the non-symmetry of the distribution. As we do not have quantitative information on the distribution for the recruitment stretch λr, we allow only a change of the interval λi and choose to fix the ratio of a1 and a2 in the fitting procedure. For all layers, the parameters for the beta function lie in the same range and below 2 indicating that the distribution of the fibre stretch is rather broad.
The mean fibre angle α is a fitting parameter. The anisotropic model is very sensitive to α preventing us from using the averaged data by Schriefl et al. , which combines the information of 11 donors. The determination of the mean fibre angle of the individual donor would be necessary. Experimental data for the fibre dispersion of thoracic aortic tissues, from which the parameters b1 and b2 are derived, are reported in . The authors fit the data with a rotationally symmetric function. We assume different concentration parameters for the in plane and the out of plane distribution as described in equation (2.12). The parameters of determination for intima, media and adventitia are R2 = 0.987, 0.951, 0.989 for these preceding fits.
The majority of the data used for fitting the model are based on experiments conducted on the thoracic aorta of elderly people. The mechanical properties, especially of the intima, change notably with age. When more experimental data are available, the model can help us to understand the changes in the microstructure of the tissue with age.
3.4. Model sensitivity
We perform the sensitivity analysis with the parameter values of the adventitia in table 1 as the reference point in the parameter space and plot the results in figure 5. The overall behaviour also applies to the media and the intima (we omit the corresponding plots). For both the circumferential and the axial direction, see figure 5a,c, the proposed model is most sensitive to α and λ1 over the whole range of λu. The sensitivity for the mean fibre angle α is negative for the circumferential direction and positive for the axial direction. This is explained by the increasing stiffness of the tissue in the axial direction when increasing the mean fibre angle with a simultaneous decrease of stiffness in the circumferential direction. The relatively high negative sensitivity owing to λ1 emphasizes the importance of the fibre contribution to the overall behaviour. Figure 5b,d zooms into the majority of the sensitivity curves in the circumferential and axial directions, respectively. The sensitivity to the shear modulus of the matrix material is intuitively maximum at the beginning and it decreases with increasing λu, whereas the opposite is true for the fibre-related parameters. It is worth noting that the sensitivity curves for μf and n match exactly. This means that these parameters are linearly dependent and at least one of them needs to be determined without fitting to create an identifiable system.
3.5. Fibre families in common iliac arteries
The data of iliac arteries, as documented by Schriefl et al. , suggest that the media comprises only one family of fibres. Here, we make a short excursus to the fit of the fibre distribution in the iliac media to show that the reported distribution can be caused by two families of fibres as well as by one family of fibres.
The experimental data and the corresponding fit of ρ are shown in figure 6; the mean fibre angle is prescribed as α = 7.5°, and the concentration parameters b1 and b2 are fitted. From the experimental data, one may conclude that there is only one family of fibres in the circumferential direction. Nonetheless, with a mean fibre angle of α = 7.5° (b1 = 29.0, b2 = 3.2), we obtain an equally good fit as for α = 0° (b1 = 29.3, b2 = 2.7) with the coefficient of determination R2 = 0.931 and R2 = 0.932, respectively.
An explanation is illustrated in figure 7 where the orientation distribution of both fibre families (ρ1 and ρ2) are plotted separately in (a) and as the sum (ρ = ρ1 + ρ2) in (b). Even though the mean fibre direction for each fibre family deviates from α = 0, the addition of both distributions gives a peak at α = 0. Figure 7c shows the coefficient of determination R2 in terms of the prescribed mean fibre angle α. The goodness of fit is equally good for angles from α = 0° to about α = 10°. Only for larger angles, the goodness of fit decreases. The results show that there might be one or two fibre families in the media of an iliac artery, and that motivates further studies of the microstructure of the iliac media.
3.6. Collagen fibre recruitment
Here, we aim to show how the collagen fibres are recruited for different deformation modes using the fitted material parameters of the adventitia as summarized in table 1. Figure 8a illustrates the ODF ρ of fibres over a unit sphere. Figure 8b–d shows the relative fibre recruitment or cumulative recruitment distribution function (CRDF), i.e. the percentage of the fibres that are recruited in each direction for an equi-biaxial, a uniaxial and a simple shear test, respectively.
Note that the absolute number of fibres recruited in a particular direction is determined through the number of fibre bundles n, the ODF ρ of fibres and CRDF. We choose stretches of λe = 1.2 and λu = 1.4 for the equi-biaxial and the uniaxial tensile test, respectively. The amount of shear γ in the simple shear test is chosen to be 0.2. For the equi-biaxial deformation, the fibres in the plane of the deformation are recruited evenly (figure 8b). For uniaxial extension in the x-direction, the fibres are recruited and aligned close to that direction (figure 8c). Owing to the incompressibility condition, we have shrinkage in the y- and z-directions with no fibre recruitment. The higher the stretch of the tissue, the smaller is the transition zone from no recruitment to complete recruitment. The contours of CRDF are circles, because the two principal stretches in the y- and z-directions are equal.
Figure 8d shows the CRDF for simple shear with shear in the x-direction. The fibre recruitment is concentrated around the Lagrangian principal strain direction. For the simple shear case, the contours do not reside in circles because simple shear, in general, possesses three different principal stretches.
We have proposed a structural and multiscale constitutive model for arteries considering fibre recruitment and three-dimensional fibre distribution. Derived from the waviness of the fibres, we have defined a recruitment stretch for the collagen fibres which exhibit a linear force–stretch relationship. For the recruitment stretch a PDF accounts for the different degrees of waviness of the collagen fibres within a fibre bundle. A bivariate von Mises distribution defines how the fibre bundles are oriented in space. These assumptions are in accordance with the morphological findings documented in the literature. We have highlighted that the variations in the properties of arteries may easily be tuned by adjusting one or more of the above-described properties.
We need more experimental data of the fibre crimping and thus the distribution of recruitment stretches of the fibres. In this study, we have treated the parameters describing the recruitment distribution as phenomenological. The different types of collagen in the layers of the arterial wall could be considered by using different Young's moduli for the different collagen types. Detailed experimental data of the distribution of each collagen type and the stiffness of the collagen fibres of the different types would be ideal for the inclusion into the proposed constitutive model. Regarding the structure of the collagen network, cross-links between fibrils play a role at the fibril-to-fibre level . This is not considered in this study and leaves room for more detailed structural models capturing the mechanical response of arteries.
In summary, the values for all parameters are within a physically explainable range. The lack of experimental results, however, prevent a more rigorous evaluation. The proposed model motivates further systematic investigations of the structural properties of soft biological tissues, in general, and of arterial tissues in particular.
This work was partly supported by the European Commission under the Seventh Framework Programme, grant agreement number no. 248782. We thank David Pierce for his support.
↵1 Version R2010b, by The MathWorks Inc., Natick, MA, USA.
- Received February 9, 2015.
- Accepted February 20, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.