Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Biomech Eng. Author manuscript; available in PMC 2011 June 6.
Published in final edited form as:
PMCID: PMC3108498

A Microstructurally-Driven Model for Pulmonary Artery Tissue


A new constitutive model for elastic, proximal pulmonary artery tissue is presented here, called the Total Crimped Fiber Model. This model is based on the material and micro-structural properties of the two main, passive, load-bearing components of the artery wall, elastin and collagen. Elastin matrix proteins are modeled with an orthotropic neo-Hookean material. High stretch behavior is governed by an orthotropic crimped fiber material, modeled as a planar sinusoidal linear elastic beam, which represents collagen fiber deformations. Collagen-dependent artery orthotropy is defined by a structure tensor representing the effective orientation distribution of collagen fiber bundles. Therefore, every parameter of the total crimped fiber model is correlated with either a physiologic structure or geometry or is a mechanically-measured material property of the composite tissue. Further, by incorporating elastin orthotropy, this model better represents the mechanics of arterial tissue deformation. These advancements result in a micro-structural total crimped fiber model of pulmonary artery tissue mechanics which demonstrates good quality of fit and flexibility for modeling varied mechanical behaviors encountered in disease states.

1 Introduction

It is important to utilize a suitable constitutive model when using the tool of structural modeling to further the understanding of pulmonary arterial hypertension (PAH). Because PAH has many effects on passive vascular wall properties and composition [19], these effects should be captured within specific parameters of the constitutive model. The extracellular matrix protein elastin determines the low-stretch stiffness of the artery, and thus the pulmonary vascular input impedance [1]. Collagen in the vessel wall limits distension [2] and also contributes to vascular stiffness. The arrangement of the fibers influences the anisotropy and stiffness of the vessel wall. Changes due to PAH in elastin matrix stiffness [3] and collagen arrangement [4] should be captured in the model. Because pulmonary arteries remodel with PAH, we seek to develop a constitutive model that captures the known changes in both pulmonary arterial elastin and collagen independently.

A constitutive model suitable for the artery wall under PAH and normal conditions should reproduce two characteristic features of vascular tissues. These features are the J-shape stress-stretch curve, observed under the application of tensile loads, and the anisotropic nature, observed when loading in different directions [2, 5, 6]. The J-shape is due to the delayed engagement of the so-called collagen fiber bundles (CFB) [79]. CFB, common structures in fibrous tissue, such as tendon, artery, and skin, in the artery wall are tortuous [8] and lightly cross-linked [2]. This feature has been modeled in the past using phenomenological models of the Fung-type [10] as well as rubber-like limited extension flexible chain models [11] successfully [12, 13], but more recently the gradual recruitment and engagement of tortuous CFB with increased stretch has been modeled [1417]. CFB have also been modeled using a crimped elastic fiber, which produces low, but non-zero forces at low stretch, but develop into linear elastic behavior at high stretches [1821]. The anisotropy of the artery tissue is due to the underlying microstructural arrangement of the CFB and elastin [2226]. The anisotropy of the overall artery wall was characterized by Patel and Fry as being orthotropic [5]. Recently, it has been noted that arterial elastin, as well as CFB, has an orthotropic nature [26, 27] yet few models have taken this into account [23]. Other models have used only orthotropic structure in the collagen portion [14, 28]. Recently, a model has been developed with a distributed collagen fiber orientation, so as to more realistically capture the true effect of out-of-plane fibers [15].

This paper presents a two-part constitutive model, representing CFB and elastin. The CFB portion of the model uses a distributed CFB orientation as in [15], but also incorporates a second tunable axis for the orientation distribution. CFB are modeled as sinusoidal elastic beams which behave linear elastically, as was found by Sasaki [21]. This approach eschews the use of an engagement distribution function, which lowers model complexity. Drawing from test data on purified PA elastin-only tissues, an orthotropic elastin model was chosen [29]. The elastin is modeled using a neo-Hookean solid [27] with neo-Hookean fiber reinforcements [30]. This approach is similar to that used in Rezakhaniha and Stergiopulos [23], but uses the neo-Hookean model for the isotropic contribution instead of a modified neo-Hookean. The paper is organized as follows: Experiments on calf tissue to obtain the data are first presented in Section 2. Based on experimental observations, the model is developed for 3D finite deformation in Section 3. In Section 4, the model is applied to uniaxial finite deformation, and is compared to the experimental data. In section 5, the effects of the parameters are explored.

2 Experimental Method

Mechanical tests of uniaxial tension of artery tissue in axial and circumferential directions were conducted. Due to the paucity of test data on large pulmonary arteries, uniaxial tests were conducted. For more detail, see [3]. To collect data for analysis and data fitting, proximal pulmonary artery tissue was taken from 2-week old hypertensive and normotensive calves and dissected to obtain samples for uniaxial strip tests. Pulmonary hypertension was induced by holding the animal at a simulated high altitude for 2 weeks, causing hypoxia leading to hypertension. Strips were cut in the axial and circumferential directions (shown in Fig. 1) from the left and right branches as well as the main trunk and the portions of the arteries with branch points, localized thickening, or other anomalies were avoided. Significant efforts were made to keep the tissue slender. However the thickness of the sample dictates the width of the strip. A too-narrow strip would cause the strip to shear and twist when the grips were tightened and can be easily damaged during sample processing. The width of the strip was usually proportional to the thickness, but in some cases, the length of the tissue obtained was not sufficient enough to allow very slender strips while keeping the width and thickness proportional. The dimensions of the strips were dependent on the branch section and calf, with most having width of around 3–8mm and lengths of 20–30mm for branch sections and 40–60mm for trunk sections. Length to width ratios were in the range of 7–12. Loose connective tissue was trimmed from the strips. Mechanical tests were conducted on a material mechanical testing machine (Insight 2, MTS) with an environmental bath filled with calcium-free phosphate buffered saline solution of pH 7.4 at 37° C. The contribution of smooth muscle cell to PA mechanical behavior was evaluated by adding rho-kinase inhibitors and contractile agents to the solution, we found no change in the mechanical behaviors of PA two hours after inoculation. Typical maximum stretching force was ~3N, varying on the sample dimensions. A five-Newton load cell afforded good resolution in this force range while still capturing the high forces experienced at higher stretches. The strain rate was kept constant at 10% per second. During testing, it was observed that the effect of the clamping seen in the transverse direction decayed rapidly along the length to a range of less than 5mm. Although the test is not as accurate as using an extensometer, the results are accurate enough to use for constitutive modeling. The tissues were tested usually within 24 hours of sacrifice, up to 48 hours after sacrifice. No change in material behavior was noted for tests performed on day one through day three post sacrifice. Drug tests were performed for smooth muscle cell contraction and inhibition, with no quantifiable changes.

Figure 1
Schematic of strips taken from excised pulmonary arteries.

The maximum stretch applied to the strips was increased gradually as needed from a uniaxial stretch of 1.3 up to 1.8 for each strip to determine the stretch necessary to engage the collagen fiber bundles. To avoid tearing, significant damage or breakage, the sample was repeatedly unloaded to note the degree of damage in the force-extension curves. Peak load was much lower than the force to tear or damage the sample. Typical maximum stretches were from 1.6 to 1.8. The specimens were then tested to that stretch for 10 cycles; nine for preconditioning, with the last cycle used as the stress-stretch data for that specimen.

3 Modeling

3.1 Preliminaries

As a preliminary to the mechanics in this paper it is necessary to define the continuum mechanics used. The conventions used in this paper largely follow Holzapfel [31]. The right Cauchy-Green deformation tensor, C, and the finger tensor b are defined as

C=FTF,  b=FFT,

where F is the deformation gradient. If the strain energy density function of an incompressible material is denoted as ψ, the second Piola-Kirchhoff (PK) stress tensor is defined as

S=2[partial differential]ψ[partial differential]CpC1,

where p is the hydrostatic pressure. The first PK stress, P, and the Cauchy stress is denoted as σ, are calculated as

P=FS,  σ=J1FSFT,

where J is the determinant of the deformation gradient. It is convenient to represent the strain energy function in terms of the invariants of C. In the isotropic case,




In the incompressible case, we require J = 1, thus the strain energy depends only on the first and second invariants.

In the anisotropic case with one family of fibers there are two additional invariants pertaining to the deformation of the fibers. If the fibers are aligned in the direction of a0 then the structure tensor following Spencer [32] is defined as

A0=a0[multiply sign in circle]a0,

and the invariants describing the deformation of the fiber family are

I4=C:A0.  I5=C2:A0.

The fourth invariant has a straightforward meaning, and can be calculated as


where λa is the fiber stretch, and the fifth invariant is related to how the fibers couple to shear deformations. In this analysis, we use a unit vector a to denote a0 in the current coordinates. Thus, a is calculated as


In the case of two fiber families, there are five additional invariants of the deformation tensor. If the fibers are aligned in the directions of a0 and g0, the structure tensors characterizing these fiber families are

A0=a0[multiply sign in circle]a0,  G0=g0[multiply sign in circle]g0.

In addition to the invariants from the transversely isotropic case, I4 and I5, there are four additional invariants. The first two are calculated similar to I4 and I5 as,

I6=C:G0,  I7=C2:G0.

The eighth and ninth invariants use both structure tensors from the fiber families, and are calculated as

I8=tr(CA0G0),  I9=tr(A0G0).

If a0 and g0 are orthogonal to one another, the eighth invariant is identically zero, and thus does not enter in to the calculations. Following Holzapfel [28], it is assumed that artery tissue is incompressible, and does not depend on the invariants involving C2. The latter is made for a number of reasons: to reduce the number of material constants and also to reduce the complexity of the model. Thus, we consider invariants I1, I4 and I6 with the constraint that J = 1.

If we consider the extension of a single fiber of length LF with cross-section area AF, deformed by observed stretch λF, the energy, ψF, is given by


where PF is the nominal stress in a fiber. With the force in the fiber, FF, equal to the nominal stress times the area, AF, the energy can be written in terms of the force generated by the fiber as


It is important to remember that when using this energy, it is not a strain energy density, but the total energy of the fiber.

3.2 Modeling an Individual Collagen Fiber Bundle

In order to model an isolated collagen fiber bundle, it is assumed that the collagen fiber bundles have a wavy configuration in the artery wall. The planar case is considered here, as was adopted previously [20, 33, 34]. Comninou and Yannas [18] developed a constitutive model for tendons consisting of sinusoidal collagen fiber bundles, termed as crimped fibers. In Comninou and Yannas [18], only small stretches and small-amplitude fiber crimp were considered. In this section of the paper, the crimped fiber model is extended to nonlinear stretch behavior. The crimped fiber model simplifies a collagen fiber bundle into a planar sinusoidal elastic beam. Because the material deformation of the collagen fibers is much less than the observed apparent deformation, the force can be calculated using a small deformation material formulation. The small deformation case holds for collagen material deformations up to 10%. The fiber has a given elastic modulus, E, cross-sectional area, A, second moment of inertia, I, period 4l0=2π/b and amplitude a as shown in Fig. 2. The parameter b is the spatial frequency of the sinusoid. The undeformed and the deformed contour lengths of the fiber are denoted as L0C, and LC, respectively. The projected length of the undeformed fiber is L0 and the deformed is L. It is important to understand the difference between the projected length and the contour length, as the average stretch internal to the fiber, defined as λc=Lc/Lc0, gives rise to the force, and the apparent stretch, defined as λF = L/L0, relates to the overall material deformation.

Figure 2
Schematic of crimped fiber model. The solid line is the undeformed fiber configuration and the dotted line is the deformed fiber configuration.

If the undeformed beam is described by


where a is the amplitude and b is the spatial frequency of the sinusoid, the total contour length of the beam is


where E(•) is the complete elliptic integral of the second kind. Using the linearization of error O(a3b3/3) applied by Comninou and Yannas [18], under the application of the tensile load, FF, the beam is assumed to deform to a sinusoidal shape with smaller amplitude but longer wavelength, expressed as


where λA is the ratio of the current amplitude to the original amplitude, calculated as


where ω relates the axial and bending deformations of the beam and is calculated as


R can be thought of as the radius of gyration for the cross section of the beam. Therefore, as R increases, the apparent bending stiffness increases. The contour length of the deformed beam can be calculated as


The undeformed contour length and undeformed projected fiber length are then related to the deformed contour length and deformed projected fiber length by


The relationship governing the stretch in the fiber compared to the apparent stretch is


Eq. 20 can be simplified and linearized with the same order of error as before, then rewritten in terms of the apparent fiber stretch, λF, and the material fiber stretch, λc, by substituting the relationships, giving


A more detailed derivation of this relationship is shown in appendix A. This relationship can be inverted to solve for the material stretch, λc, given λF. Assuming the fiber behaves as a Hookean material, given the observed stretch, the force is calculated from the stretch as


Here, the fiber is treated as a linearly elastic beam because fiber strains are low. Although the overall material deformation can be large, the fiber material does not undergo high strain due to rigid body rotation of fiber and bending deformation, which can accommodate large deformation. Large material deformations of the fibers are not usually experienced under normal operating conditions, and as such, a small fiber strain measure can be adopted. The small fiber material strain can be seen in the Session 4.2. From Fig. 2, the six parameters, A, E, a, R, l0 and I can be further simplified to four parameters: the modulus, E, the cross-sectional area, A, the ratio between bending rigidity and extensional stiffness, R2/l02=4I/Al02, and another specifying the geometry of the fiber, [theta w/ overline]0 = atan(2ab/π). Using these new parameters and rewriting Eq. 21 give


Using this equation for the calculation of force from stretch, the behavioral dependence on the stress can be observed. The change in fiber behavior due to varying R/l0 while holding [theta w/ overline]0 constant is shown in Fig. 3(a). At low R/l0 values (R/l0 <0.1) there is very little stiffness in the low-stretch region. With increasing R/l0, the fiber bears load earlier, but at high stretches retains the same fully-developed apparent stiffness. This parameter thus determines the low-stretch stiffness and the shape of the transition from low-stiffness to fully-developed stiffness.

Figure 3Figure 3
The behavior of a single fiber under uniaxial extension is plotted. The nominal fiber stress, P is normalized by the Young’s modulus to give the reduced fiber stress. A) As the bending stiffness, via R/l0, is increased, the fiber behavior goes ...

In Fig. 3(b), [theta w/ overline]0 is varied while holding R/l0 constant for a single, isolated fiber under uniaxial tension. The largest effect of increasing this parameter is to increase the observed stretch at which the fiber begins to bear load. In addition to increasing the engagement stretch, [theta w/ overline]0 changes the behavior of the engagement. For example, when [theta w/ overline]0 is small, a sharp transition of slope can be seen. The transition broadens with increasing [theta w/ overline]0. In addition, because [theta w/ overline]0 also increases the total length of the fiber it has the added effect of decreasing the fully-developed apparent modulus even though the material modulus has not changed.

This fiber model captures the transition from low stiffness to high stiffness using two observable parameters. The parameters can be obtained from histology and applied to the fiber model or can be fit to data taken from tissues using a least squares method. By using histological measurements of the fiber shape, the tortuosity can be measured, and applied to the sinusoid in order to obtain [theta w/ overline]0. This method is to be developed in the future.

3.3 Orthotropic Crimped Fiber Model

Because artery tissue can be thought of as locally orthotropic, it is appropriate to utilize an orthotropic constitutive formulation. To this end, a three-dimensional structure tensor, with two orthogonal directions and an isotropic component is used to characterize the material. As with the distributed collagen fiber orientation method of modeling used previously [15], a three-dimensional structure tensor is utilized to characterize the orientational distribution of the collagen fiber bundles. Unlike the method employed by Gasser et al [15], the structure tensor used here has two orthogonal groups of fibers, forming the major axes of an ellipsoid with tunable shape in all three major directions. The shape of the ellipsoid represents the angular distribution of collagen fibers in the artery wall. Collagen fibers in the artery wall have a distribution of orientations. The mean direction of orientation follows the circumferential direction. Similar to the symmetric fiber families about the circumferential direction used in [28], which assumes that the fibers exist in the circumferential-longitudinal plane, a second direction is needed to define the plane in which the fibers lie and the local directions governing the orthotropy of the tissue. We define the local coordinate system through the vectors describing the circumferential and longitudinal directions. In addition to the vector describing the circumferential direction, the second vector is necessary to form the plane. In this paper, the second vector is chosen to be the longitudinal direction. Thus, the orthogonal vectors a0 and g0 are aligned with the circumferential and longitudinal directions respectively, defining the plane in which the fibers generally lie. The structure tensor then takes the form

H0=κ+γ13I+(1κ)a0[multiply sign in circle]a0+(1γ)g0[multiply sign in circle]g0,

where κ and γ are the structure parameters which affect the shape of the tensor and must be between zero and one, and also satisfy κ + γ ≥ 1, which ensures that the dimensions of the structure tensor are positive. A negative eigenvalue in the structure tensor would cause the model to lose convexity and violate the global minimum of strain energy at no deformation. This effect is detailed in Appendix B. This structure tensor can be visualized as an ellipsoid, with trace of 1, as in Fig. 4. When the structure parameters are both unity, the coefficients for the second and third terms in Eq. 24 vanish and the structure tensor represents a sphere, and the behavior of the material is isotropic. The structure tensor’s shape reflects the fiber distribution and gives rise to the model’s anisotropy. The two parameters in the structure tensor currently are fit to the test data, but we are currently working toward histological image and structure parameter correlation. Typical parameter values for κ and γ are around 0.8–0.9, which gives an oblate ellipsoid with thickness around 0.2. The structure tensor in the current configuration is denoted as


As with the model introduced by Gasser et al. [15] and Spencer’s work [35], the stretch experienced by a fiber is related to the structure tensor by λF2=H0:C. Using Eqs. 6 and 10a, this can be written in terms of the invariants,

Figure 4
The visualized ellipsoidal structure tensor. This is generated by observing how the structure tensor transforms a unit vector. A longer dimension indicates a higher concentration of fibers in that direction. Note here that a0 and g0 are not aligned with ...

The strain energy density function is given by


where K is the density of the fibers in the composite. In general, for a material dependent on the first, fourth and sixth invariants, the second PK stress is given by

SCF=2[partial differential]ψCF[partial differential]I1IpC1+2[partial differential]ψCF[partial differential]I4a0[multiply sign in circle]a0+2[partial differential]ψCF[partial differential]I6g0[multiply sign in circle]g0.

The relevant terms to compute are the derivatives of the strain energy function with respect to the invariants,

[partial differential]ψCF[partial differential]I1=12KFF(λF)λF1(κ+γ13)[partial differential]ψCF[partial differential]I4=12KFF(λF)λF1(1κ)[partial differential]ψCF[partial differential]I6=12KFF(λF)λF1(1γ).

Using Eqs. 28 and 29, the second PK stress becomes

SCF=KFF(λF)λF[(κ+γ13)I+(1κ)a0[multiply sign in circle]a0+(1γ)g0[multiply sign in circle]g0]pC1.

In the isotropic case, κ = γ = 1, the strain energy function should reduce to an isotropic model and the second PK stress is


The second PK stress and the Cauchy stress from Eq. 30 can be rewritten using Eq. 24 respectively as



3.4 The orthotropy of arterial elastin

From elastin test data taken from purified pulmonary artery samples [29], as well as thoracic aorta sections [36], it was seen that the behavior of elastin differed in different directions. Therefore, anisotropic material model is used for elastin. For the model on elastin and non-protein matrix, we assume it consists of two orthogonal families of neo-Hookean fibers oriented with the axial and circumferential directions are embedded in a neo-Hookean matrix. As in deBotton and Socolsky’s model for a neo-Hookean fiber [30], we use the strain energy function for the fiber reinforcement, ψEla of


where I4 is the fourth invariant associated with the circumferential direction, which also coincides with the a fiber direction. The effective modulus of the fiber is characterized by μa. We take the derivative with respect to I4,

[partial differential]ψEla[partial differential]I4=μa2(1I432),

which is used in the stress formulation. Unlike the model presented by Rezakhaniha and Stergiopulos, which used a transversely isotropic material [23], we choose an orthotropic representation. Together with the isotropic neo-Hookean matrix of shear modulus μ and an additional fiber family coincident with the axial direction with shear modulus μg, the total strain energy for the elastin component is given as


The second PK stress and the Cauchy stress are

SEl=μI+(μaμ)(1I432)a0[multiply sign in circle]a0+(μgμ)(1I632)g0[multiply sign in circle]g0pC1.

σEl=μb+(μaμ)(I4I4½)a0[multiply sign in circle]a0+(μgμ)(I6I6½)g[multiply sign in circle]gpI.

This model can easily be reduced to the isotropic neo-Hookean case by setting μa and μg to μ. The ratio between the circumferential and the longitudinal moduli can be large, as it has been found that thoracic aorta has circumferential modulus twice that of the longitudinal direction [36]. Extreme anisotropy, a ratio of 10 or more, is not expected, and typical anisotropy ratios range from 1 to 5.

3.5 Complete model incorporating volume fractions

The model strain energy is the total of the strain energies representing the collagen fiber bundles and the elastin network weighted by its volume fraction,


Since the collagen fiber bundle model already contains a density (K), Eq. 37a does not have a volume fraction for collagen fiber bundle energy. In addition, it is possible to lump the volume fraction fEl with elastin shear moduli. This linear superposition of strain energy allows us to calculate the stresses as


which can be rewritten as

S=[μI+(μaμ)(1I432)a0[multiply sign in circle]a0+(μgμ)(1I632)g0[multiply sign in circle]g0]+KFF(λF)λFH0pC1.

Because A and K together determine the total area of the collagen fibers per unit material area, it is advantageous to lump the two together, KA. In total, there are nine parameters to consider here: three moduli concerning the elastin protein network, isotropic shear modulus, μ, circumferential shear modulus, μa, and longitudinal shear modulus, μg; three parameters for the collagen fiber bundle; intrinsic collagen fiber Young’s modulus, E, fiber shape, [theta w/ overline]0, and normalized radius of gyration, R/l0; one pertaining to collagen fibers per unit material area, KA; and two for the orthotropic structure tensor, circumferential, κ and longitudinal, γ. For the following sections, the intrinsic material modulus of the collagen fiber, E, was chosen to be in the range of previous work [14, 21, 37, 38], with a value of 10 GPa. The model from here on will be called the total crimped fiber (TCF) model.

4 Results

In the following, we evaluate the TCF model behavior by investigating the stress-strain response of artery tissue under uniaxial loading conditions. A more thorough evaluation of the model is to investigate its performance under different loading conditions, such as uniaxial and biaxial. This is currently attempted by the authors and will be reported in the future. First, the stress-stretch response of the model under uniaxial loading is developed. Second, we conduct parametric studies to observe the effects of structural parameters on the model behavior. Third, we compare the model behavior with experimental results. Finally, we use the model to explore parameters used previously in judging artery behavior.

4.1 Model responses under uniaxial deformation

Material point simulations were performed using the model with varying parameters. Applying the load in the 1-direction, S11, one can solve for the uniaxial stress by applying the appropriate boundary conditions of S22=S33=0 and the incompressibility constraint. Taking the full orthotropic model, with the stress applied in the 1-direction,




Using S33=0, we can solve for the Lagrange multiplier associated with the incompressibility constraint, p, as


Combining this with S22, we obtain the expression


This can be solved numerically for λ2 given λ. Once λ2 is obtained, λF and the uniaxial stress can be calculated. For uniaxial tests in the 2-direction where S22 is controlled, the process is similar, but with S11=0. The uniaxial deformation case was used to fit the experimental data in the following sections.

4.2 Comparison with Experiments

The model was fit to the uniaxial data taken from the excised calf tissues to evaluate the model both qualitatively and quantitatively. Some data fits are shown in Fig. 5. The markers are experimental data points, while the lines are the model fits. As seen here, the fit quality is very good in different situations of material behaviors. In Fig. 5A, it is seen that the material behaves nearly transversely isotropically in the low-stretch regions, and grow apart with increased stretch. This is a function of the anisotropy of the collagen fiber portion of the model. In Fig. 5B, the lower-stretch portion of the behavior have different stiffnesses, but this is consistent with the engagement of the collagen fibers. To fit these data, the anisotropic emphasis is placed on the collagen portion. In Fig. 5C, we see a nearly transversely-isotropic response, and this is reflected in the parameters in table 1. Fig. 5D. shows an interesting behavior in which the low-stretch behavior is highly anisotropic while the collagen engagement stretch does not change. The transition of collagen engagement is captured well by the model. The high-stretch portions of the curves are also fit well. The material parameters used to fit the data are presented in Table 1. While the parameters corresponding to the neo-Hookean portion vary considerably, it is important to note that the low-stretch behavior does vary from tissue to tissue. Though the isotropic neo-Hookean portion of the model is highly variable, it is consistently low. This variability is, in part, due to the use of two uniaxial tests. Better results might be obtained when using data from biaxial tests. The R2 values presented are calculated as


where y is the stress for the circumferential and longitudinal direction uniaxial tests. It is seen from the high R2 in table 1 that the model fits the data well. It captures the anisotropy of both the low-stretch portion and the crimped fiber portion.

Figure 5Figure 5Figure 5Figure 5Figure 5
The data from four sets of uniaxial tests and their corresponding model fits. Squares and triangles denote circumferential and axial data respectively, while Solid and dashed lines denote circumferential and axial data fits. Here it is seen that the model’s ...
Table 1
Fit parameters and normalized error for selected tissue sections. These correspond to figures 5(A–D) and 8(A–D).

Fig. 5E shows the fiber material stretch as a function of the material stretch for the case of circumferential test in Fig. 5B. The fiber material stretch is calculated from Eq. 20. It is seen that as the material stretch reaches ~1.6, the fiber material stretch is only ~1.009. As discussed before, this is due to the fiber rigid body rotation and fiber bending can accommodate a large amount of stretch before the fiber itself being stretched largely.

5 Discussion

5.1 Parametric Study

The model has nine independent parameters. The orthotropic neo-Hookean portion of the model is quite straight-forward in its behavior. The parameter μ defines the out-of-plane stiffness in the low-stretch region. This parameter comes into play when the material is deformed biaxially. It determines the coupling between the stresses in the a, or circumferential, and g, or longitudinal, directions through the incompressibility constraint and the prescribed boundary conditions. A larger μ will cause a larger stress to be transmitted to the orthogonal direction for a given deformation state, and vice versa. The parameters μa and μg drive the low-stretch stiffness under uniaxial deformation in the circumferential and axial directions respectively. Changes in these parameters indicate a change in the elastin network. Any increase here indicates an increase in low-stretch stiffness, which is an observed change in pulmonary hypertension [3]. The strongest indicator of large artery stiffness is μa, as it determines the circumferential stiffness, and therefore affects the volumetric capacitance.

The crimped fiber portion of the model has parameters that pertain to the fibers themselves, and parameters that pertain to the orientation of the fibers. The effects of the collagen portion of this model can be elucidated by how they relate to the material microstructure. The collagen fiber number fraction K, along with the Young’s modulus of the fibers, E, and single-fiber area, A, define the stiffness of the crimped fiber portion of the model. As we have shown earlier, it is conceptually simple to fix the modulus and lump the number density and fiber area into an effective area fraction, KA. The structure tensor characterizes the degree of anisotropy. The shape of the sinusoidal beam, [theta w/ overline]0, drives the engagement point of the fibers. It relates to the tortuosity in the collagen fiber bundle, and thus drives the engagement point of the fibers, relating to the crimp measured by Hurschler et al[39]. An increase in this parameter will cause an increase in the engagement point and a subsequent decrease in the fully-developed stiffness due to the additional total length of the fiber. This would indicate that elastin would continue to dominate at a higher stretch than with a lower value of [theta w/ overline]0. The bending stiffness parameter, R/l0, determines the breadth of the transition between elastin-dominated low-stretches and the collagen-dominated high-stretch regimes. An increase here would indicate either: a fiber with larger bending stiffness, a change in the fiber properties; or a broader distribution of the engagement stretch of the collagen fibers.

In the following section, the effects of the crimped fiber orthotropic model parameters are studied through a parametric study. The parameters can be tuned to achieve a certain behavior. Besides having high flexibility in modeling anisotropy and engagement, the model produces a basic shape of the stress-stretch curves that is consistent with the J-shape common to artery behavior which simplifies the use and understanding of the model.

Fig. 6(A–C) show the stress-stretch curves for varying one parameter and holding the others constant. In these figures, λ1 and λ2 correspond to the axial and circumferential stretch, respectively, while λ3 is the radial stretch. In calculating the material response, the uniaxial stretch, λ1 or λ2, is the independent parameter, and σ1, λ2 and λ3 or σ2, λ1 and λ3 are the dependent variables, as in Eq. 39. In Fig. 6(A–C), the uniaxial material response is shown for the orthotropic crimped-fiber model only, showing the effect of changing the model parameters pertaining to the fibers and their orientations. This was performed by setting the three shear moduli of the orthotropic neo-Hookean components to zero. For each set of parameters, the resultant uniaxial behaviors along the two fiber directions are shown.

Figure 6Figure 6Figure 6
A) A parametric study varying γ, with κ = 0.8. It is seen that when κ = γ, (shown in circles), the behavior is transversely isotropic; the uniaxial stress-stretch curves lay on top of one another. As γ is increased ...

In Fig. 6(A) the change in material behavior due to changing γ while holding all other parameters constant causes an observed change in the anisotropy of the collagen component. A higher value for γ corresponds to less fiber alignment in the longitudinal direction, and thus the stiffness and fiber engagement stretch in the longitudinal direction is directly affected by this change. The same is true for κ and the circumferential direction. If both parameters decrease in the same proportion, the structure tensor flattens, and causes a decrease in the coupling between the transverse direction and the in-plane directions. A highly anisotropic material can be specified with the given parameters, but the fully developed stiffness will differ with direction.

The change in material behavior due to a change in [theta w/ overline]0, as in Fig. 6(B), is observed in both the ultimate stiffness and engagement stretch. Increasing [theta w/ overline]0 directly causes an increase in the engagement stretch. In addition to increasing the engagement stretch, the ultimate stiffness is decreased. This is, as stated earlier, due to the increased contour length of the beam with constant end-to-end distance. It may appear that the degree of anisotropy is increased with increasing [theta w/ overline]0, but the perceived increase in degree of anisotropy is due only to the increased stretch at which engagement occurs.

The effect of R/l0 on behavior is explored in Fig. 6(C). At low values for R/l0, there is very little stiffness at low stretches, which rapidly increases and becomes linear quickly above the engagement stretch. At higher values, it is seen that there is some low-stretch stiffness and the transition to the fully-developed stiffness is more gradual.

While the parameters are not orthogonal, that is, a change in one parameter may change more than one observable characteristic, the parameters are physically relevant. The change of the shape of the fiber through [theta w/ overline]0 may change the fully-developed stiffness, but it does so in a physically-based way. An orthogonal set of parameters could be found in order not to change aspects of the behavior not directly associated with the parameter in order to more easily understand the behavior, but those parameters would not be physically-based, and so would not directly correlate to any microstructural characteristics. Some parameters of this model are not unique, as some parameters, such as the fiber density, K, and collagen modulus, E, would affect the model in the same way. However, it is conceivable that two or more microstructural changes can affect the same changes in the mechanical behavior under certain mechanical loading conditions (see [14] for a brief discussion on uniqueness of collagen fiber stiffness and engagement). One way to overcome this deficiency is to connect model parameters with certain features of the stress-strain curves and then to use these connections to guide material parameter identification. However, this approach could become cumbersome if the number of experiments from different tissues increases. The other approach is to use data from experiments under different loading conditions, such as biaxial experiments with different loading ratios between longitudinal and circumferential directions.

5.2 Collagen engagement point and material parameters

Using a method described by Lammers et al [3], the collagen engagement point was correlated to the parameters [theta w/ overline]0 and R/l0. Briefly, the curvature of the stress-stretch data is calculated, and the point of maximum curvature is taken to be the collagen engagement point. Note that this is not the collagen transition point used in Lammers et al, but the engagement point designates the stretch at which the majority of collagen is engaged. The full model was used to generate the curvature data, but all parameters were kept constant save for the parameter in question. Fig. 7(A) shows engagement stretch as a function of the parameter R/l0. It seen here that the bending stiffness does affect the engagement stretch, but quickly loses sensitivity as R/l0 grows larger than 0.5. The region of greatest sensitivity is around R/l0 =0.25, which is much greater than the value for the data fits. Fig. 7(B) shows the engagement stretch as a function of [theta w/ overline]0. Though not shown on this plot, the engagement stretch would reach unity when [theta w/ overline]0 reaches zero, and would asymptote to infinity as [theta w/ overline]0 approaches 90°. The ranges of the parameters for R/l0 and [theta w/ overline]0 are around 0.001 to 0.01 and 40–50° respectively. In this range, the engagement stretch is not sensitive to R/l0, but is most sensitive to [theta w/ overline]0. Therefore it can be said that shape of the fiber affects the engagement stretch of the model more than the bending stiffness.

Figure 7Figure 7
The engagement stretch as a function of (A) R/l0 and (B) [theta w/ overline]0. The parameters used to calculate these data are μ = μa = μg =5kPa, KA=12×10−4, E=10GPa, [theta w/ overline]0 =45°, κ=0.90 ...

5.3 Model application to physiological loading conditions

The goal of this research is to present a structure-driven constitutive model for mechanical behaviors of arterial tissue. The model fits in the above sessions were based on uniaxial experimental data, where two separate samples were cut from the same animal in each data set. In order to investigate the model performance under physiological conditions, we use our recent preliminary data obtained from planar biaxial tests. Here, a planar biaxial tester was used [40]. The sample was cut from the MPA with a size of ~25mm (w) × 25mm (h) × 4mm (t) and was tested in the similar solution environment as for uniaxial experiments. Three experiments were conducted on the same tissue with the following circumferential vs. longitudinal (C:L ratio) stress ratios: 100:0 (circumferential uniaxial), 0:100 (longitudinal uniaxial), and 100:25. The circumferential and longitudinal experimental data were used to fit the material parameters. Fig. 8A shows very good fits (material parameters are show in Table 1). These material parameters were then used to predict the 100:25 biaxial stress experiment, which is close to physiological loading conditions. It can be seen that the model using material parameters from circumferential and longitudinal uniaxial test fits predicts reasonably well the 100:25 case. Fig. 8B shows that the model captures the low stretch region very well but starts to deviate from the experimental curves gradually as the stretch increases, indicating a future improvement of the model is necessary. However, we also note that the circumferential uniaxial experiment and 100:25 (C:L) biaxial experiments are closer to the physiological conditions, we can fit the material parameters using the curves from these two experiments. Figs. 8C and 8D show the model fit (material parameters are show in Table 1). It is clear that as we focus only on the physiological relevant experimental data, the model can fit the curves very well. Therefore, from the preliminary biaxial experimental data, the model is able to capture the biaxial deformation behavior of PA tissue under physiological conditions. This certainly needs to be further verified with more experimental data. On the other hand, we also realize that from a constitutive modeling point of view, one would expect that the model can capture the material behavior under a wider spectrum of loading condition, such as the biaxial loading condition with C:L ratios such as 100:0, 100:25, 100:50, 100:100, 50:100, 25:100, 0:100. We will also refine the model in order to better reflect the material behavior and structure. This is currently investigated by the authors and will be reported in the future.

Figure 8Figure 8Figure 8Figure 8
Fits from planar biaxial data. A) Uniaxial stress-stretch data from the planar biaxial tests are shown with corresponding model fits. B) The stress-stretch data from the 100:25 (C:L) experiment is shown with corresponding model prediction, using parameters ...

6 Conclusion

An orthotropic constitutive law, called the Total Crimped Fiber Model, was presented for artery tissues. The model incorporates orthotropic behavior for both elastin and collagen components, allowing the model to produce different orthotropic behaviors for the two major components. The elastin portion is a neo-Hookean fiber reinforced neo-Hookean solid. This reproduces the neo-Hookean behavior under uniaxial deformations, which agrees with the behavior as stated in literature. The collagen portion uses a microstructural basis for the nonlinear J-shape of the stress-stretch curve. A planar sinusoidal linear elastic beam is employed to model the tortuous collagen fiber bundles. The model fits the data well, and has less complexity than other models, due to the flexibility afforded by the sinusoidal crimped fiber model and the ellipsoidal structure tensor used to represent the orthotropic behavior of the material. The model can also elucidate certain aspects of the vascular behavior such as increased stiffness of either the collagen or elastin components with physiological relevant parameters. The progression of this work is to quantifiably correlate the model material properties with histological images and measurements from the material itself. Based on limited biaxial planar experimental data, the model shows that when using the material parameters from circumferential and longitudinal experimental data to predict the loading conditions close to physiological conditions, the model captures the overall behavior well and predicts well the initial loading portion but presents a stiffer response at higher stretch. However, when the experimental data from the circumferential uniaxial loading and 100:25 C:L ratio loading conditions, the model is able to fit both experimental data very well. Further work will be performed in order to fully determine the multiaxial performance of the model, and if necessary, fine-tune the model.


This work was partially supported by grants from the NIH (T32-HL072738, SCCOR-HL081506, K24-HL084923). We are also grateful for the discussions with Prof. Sacks on fiber tortuosity.

Appendix A: Derivation of contour length of the undeformed and deformed sinusoidal beam

To calculate the contour length of the beam, we find the incremental arc length ds for a given differential length dx and subsequent differential change dy, with Eq. 13 as the function,


Integrating along the length, we obtain


Integrating only to the first quarter wavelength is necessary, as the quarter wavelengths shapes are similar, we obtain


For the deformed beam whose shape is described by Eq. 15, the arc length s is calculated as


and for the first quarter wavelength is


Appendix B

In order to assess the consequences of a negative value in the structure tensor, we assume a transversely isotropic material (κ = γ). With κ = γ, and an applied stretch, λ, in the transverse direction, the deformation tensor and structure tensor are


As in equation 26, the fiber stretch is calculated as


The term within the radical is calculated as


With increasing applied stretch, λ2 increases faster than 1/λ decreases, and as such, with 2κ13<0, the term within the radical becomes negative, which leads to an imaginary number for the fiber stretch. The same effect occurs if any other dimension of the ellipsoid is negative as well. An ellipsoid dimension of 0 is permitted, but negative dimensions will give imaginary results.


1. Hunter KS, Lee PF, Lanning CJ, Ivy DD, Kirby KS, Claussen LR, Chan KC, Shandas R. Pulmonary vascular input impedance is a combined measure of pulmonary vascular resistance and stiffness and predicts clinical outcomes better than pulmonary vascular resistance alone in pediatric patients with pulmonary hypertension. American Heart Journal. 2008;155(1):166–174. [PMC free article] [PubMed]
2. Humphrey JD. Cardiovascular solid mechanics : cells, tissues, and organs. New York: Springer; 2002.
3. Lammers SR, Kao PH, Qi HJ, Hunter K, Lanning C, Albietz J, Hofmeister S, Mecham R, Stenmark KR, Shandas R. Changes in the structure-function relationship of elastin and its impact on the proximal pulmonary arterial mechanics of hypertensive calves. Am. J. Physiol.-Heart Circul. Physiol. 2008;295(4):H1451–H1459. [PubMed]
4. Jacob MP. Extracellular matrix remodeling and matrix metalloproteinases in the vascular wall during aging and in pathological conditions. Biomed. Pharmacother. 2003;57(5–6):195–202. [PubMed]
5. Patel DJ, Fry DL. Elastic Symmetry of Arterial Segments in Dogs. Circulation Research. 1969;24(1) pp. 1-&. [PubMed]
6. Doyle JM, Dobrin PB. Finite Deformation Analysis of Relaxed and Contracted Dog Carotid-Artery. Microvascular Research. 1971;3(4) pp. 400-&. [PubMed]
7. Roach MR, Burton AC. The Reason for the Shape of the Distensibility Curves of Arteries. Canadian Journal of Biochemistry and Physiology. 1957;35(8):681–690. [PubMed]
9. Elbischger PJ, Cacho R, Bischof H, Holzapfel GA. Modeling and characterizing collagen fiber bundles. IEEE; Piscataway, NJ, USA. 2006. p. 1280.
10. Fung YC, Fronek K, Patitucci P. Pseudoelasticity Of Arteries And The Choice Of Its Mathematical Expression. Am. J. Physiol. 1979;237(5):H620–H631. [PubMed]
11. Bischoff JE, Arruda EM, Grosh K. A microstructurally based orthotropic hyperelastic constitutive law. Journal of Applied Mechanics-Transactions of the Asme. 2002;69(5):570–579.
12. Zhang YH, Dunn ML, Hunter KS, Lanning C, Ivy DD, Claussen L, Chen SJ, Shandas R. Application of a microstructural constitutive model of the pulmonary artery to patient-specific studies: Validation and effect of orthotropy. Journal of Biomechanical Engineering-Transactions of the Asme. 2007;129(2):193–201. [PMC free article] [PubMed]
13. Zhang YH, Dunn ML, Drexler ES, McCowan CN, Slifka AJ, Ivy DD, Shandas R. A microstructural hyperelastic model of pulmonary arteries under normo- and hypertensive conditions. Annals of Biomedical Engineering. 2005;33(8):1042–1052. [PubMed]
14. Zulliger MA, Fridez P, Hayashi K, Stergiopulos N. A strain energy function for arteries accounting for wall composition and structure. Journal of Biomechanics. 2004;37(7):989–1000. [PubMed]
15. Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of the Royal Society Interface. 2006;3(6):15–35. [PMC free article] [PubMed]
16. Cacho F, Elbischger PJ, Rodriguez JF, Doblare M, Holzapfel GA. A constitutive model for fibrous tissues considering collagen fiber crimp. International Journal of Non-Linear Mechanics. 2007;42(2):391–402.
17. Sacks MS. Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar-collagenous tissues. Journal of Biomechanical Engineering-Transactions of the Asme. 2003;125(2):280–287. [PubMed]
18. Comninou M, Yannas IV. Dependence of Stress-Strain Nonlinearity of Connective Tissues on Geometry of Collagen-Fibers. Journal of Biomechanics. 1976;9(7):427–433. [PubMed]
19. Sasaki N, Odajima S. Elongation mechanism of collagen fibrils and force-strain relations of tendon at each level of structural hierarchy. Journal of Biomechanics. 1996;29(9):1131–1136. [PubMed]
20. Buckley CP, Lloyd DW, Konopasek M. On the Deformation of Slender Filaments with Planar Crimp: Theory, Numerical Solution and Applications to Tendon Collagen and Textile Materials. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences. 1980;372(1748):33–64.
21. Sasaki N, Odajima S. Stress-strain curve and Young's modulus of a collagen molecule as determined by the X-ray diffraction technique. Journal of Biomechanics. 1996;29(5):655–658. [PubMed]
22. Wolinsky H, Glagov S. Structural Basis for Static Mechanical Properties of Aortic Media. Circulation Research. 1964;14(5) pp. 400-&. [PubMed]
23. Rezakhaniha R, Stergiopulos N. A structural model of the venous wall considering elastin anisotropy. Journal of Biomechanical Engineering-Transactions of the Asme. 2008;130(3):11. [PubMed]
24. Sherebrin MH. MECHANICAL ANISOTROPY OF PURIFIED ELASTIN FROM THE THORACIC AORTA OF DOG AND SHEEP. Can. J. Physiol. Pharmacol. 1983;61(6):539–545. [PubMed]
25. Lillie MA, Gosline JM. Unusual swelling of elastin. Biopolymers. 2002;64(3):115–126. [PubMed]
26. Lillie MA, David GJ, Gosline JM. Mechanical role of elastin-associated microfibrils in pig aortic elastic tissue. Connective Tissue Research. 1998;37(1–2):121–141. [PubMed]
27. Gundiah N, Ratcliffe MB, Pruitt LA. Determination of strain energy function for arterial elastin: Experiments using histology and mechanical tests. Journal of Biomechanics. 2007;40(3):586–594. [PubMed]
28. Holzapfel GA, Gasser TC, Ogden RW. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity. 2000;61(1–3):1–48.
29. "Unpublished results.".
30. deBotton G, Hariton I, Socolsky EA. Neo-Hookean fiber-reinforced composites in finite elasticity. Journal of the Mechanics and Physics of Solids. 2006;54(3):533–559.
31. Holzapfel GA. Nonlinear solid mechanics : a continuum approach for engineering. Chichester ; New York: Wiley; 2000.
32. Spencer AJM. Continuum Physics volume 1: Mathematics. London, UK: Academic; 1971. Theory of invariants; p. 239.
33. Basu AJ, Lardner TJ. Deformation Of A Planar Sinusoidal Elastic Beam. Zeitschrift Fur Angewandte Mathematik Und Physik. 1985;36(3):460–474.
34. Garikipati K, Goktepe S, Miehe C. Elastica-based strain energy functions for soft biological tissue. Journal Of The Mechanics And Physics Of Solids. 2008;56(4):1693–1713.
35. Spencer AJM. Continuum Theory of the Mech of Fibre-Reinf Compos. Vienna, Austria: Springer Verlag; 1984. CONSTITUTIVE THEORY FOR STRONGLY ANISOTROPIC SOLIDS; p. 1.
36. Zou Y, Zhang YH. An Experimental and Theoretical Study on the Anisotropy of Elastin Network. Annals of Biomedical Engineering. 2009;37(8):1572–1583. [PubMed]
37. Cusack S, Miller A. Determination Of The Elastic-Constants Of Collagen By Brillouin Light-Scattering. J. Mol. Biol. 1979;135(1):39–51. [PubMed]
38. Harley R, James D, Miller A, White JW. Phonons And Elastic-Moduli Of Collagen And Muscle. Nature. 1977;267(5608):285–287. [PubMed]
39. Hurschler C, Provenzano PP, Vanderby R. Application of a probabilistic microstructural model to determine reference length and toe-to-linear region transition in fibrous connective tissue. Journal of Biomechanical Engineering-Transactions of the Asme. 2003;125(3):415–422. [PubMed]
40. Sacks MS. A method for planar biaxial mechanical testing that includes in-plane shear. Journal of Biomechanical Engineering-Transactions of the Asme. 1999;121(5):551–555. [PubMed]