## Abstract

The mechanical behavior of blood vessels is known to be viscoelastic rather than elastic. The functional role of viscoelasticity, however, has remained largely unclear. The hypothesis of this study is that viscoelasticity reduces the stresses and strains in the vessel wall, which may have a significant impact on the fatigue life of the blood vessel wall. To verify the hypothesis, the pulsatile stress in rabbit thoracic artery at physiological loading condition was investigated with a quasi-linear viscoelastic model, where the normalized stress relaxation function is assumed to be isotropic, while the stress-strain relationship is anisotropic and nonlinear. The artery was subjected to the same boundary condition, and the mechanical equilibrium equation was solved for both the viscoelastic and an elastic (which has a constant relaxation function) model. Numerical results show that, compared with purely elastic response, the viscoelastic property of arteries reduces the magnitudes and temporal variations of circumferential stress and strain. The radial wall movement is also reduced due to viscoelasticity. These findings imply that viscoelasticity may be beneficial for the fatigue life of blood vessels, which undergo millions of cyclic mechanical loadings each year of life.

- artery
- cyclic loading
- fatigue life
- viscoelastic

biological soft tissues are known to have highly nonlinear stress-strain relations, strongly anisotropic mechanical properties, and significant viscoelastic features (10). Blood vessels are typical soft tissues, and their mechanical behavior is imperative for understanding vascular processes under physiological and pathological conditions. The pseudoelastic formulation (i.e., loading and unloading follow different elastic stress-strain curves) has been widely utilized to investigate salient mechanical characteristics of blood vessels (e.g., Refs. 5, 7, 8, 10, 13). In more realistic analyses, the viscoelastic behavior of arterial walls has been studied by many researchers (4, 10, 12, 20, 21, 23, 24, 27, 29). To reduce the mathematical complexity, the theory of quasi-linear viscoelasticity (10) has been applied to arteries and other living tissues (e.g., Refs. 3, 6, 11, 15, 26). One known merit of viscoelasticity in blood vessels is to reduce the wall stress and strain during a sudden increase of mechanical loading, such as in acute hypertension (1). Although this is an interesting protective acute mechanism, the functional role of viscoelasticity in blood vessels and other biological tissues has not been fully understood.

From a mechanical point of view, the energy dissipation due to viscoelasticity may produce heat to maintain isothermal conditions in tissues, and the force-deformation hysteresis loop (or phase delay between stress and strain) caused by internal friction can filter out instantaneous changes in loading, which may prevent mechanical failure. Similarly, it has been noted that the energy storage, transmission, and dissipation in connective tissues are important for arterial performance and prevention of premature mechanical failure (21).

The residual stress and strain in unloaded arteries can homogenize the wall stress and strain at physiological condition (7, 8, 25). The vascular smooth muscle tone plays a similar role to reduce the transmural gradients of arterial wall stress and strain (1, 14, 17, 19). In the present study, we will show that the viscoelasticity of arteries can reduce the magnitudes and temporal variations of wall stress and strain at physiological loading. This result may have significant implications for fatigue life, since arteries are subjected to continuously cyclic pressure loading during the entire lifespan. Based on the model predictions, we hypothesize that the viscoelastic behavior of blood vessels may play an important role in mechanical homeostasis.

## MATERIALS AND METHODS

#### Mathematical model.

To make the problem analytically tractable, the artery is assumed to be a thick-walled tube made of homogeneous incompressible material with cylindrically orthotropic mechanical properties. The zero-stress state is represented by a sector with an opening angle Φ (Fig. 1*A*). The loaded state is approximated by an axisymmetric configuration with a uniform circular cross section (Fig. 1*B*). At zero-stress state, the inner radius *R*_{i} and outer radius *R*_{o} can be computed as (1) where χ = π/(π − Φ), and *L*_{i} and *L*_{o} denote inner and outer circumferences of the open sector, respectively.

Neglecting flow-induced shear deformation (the shear stress and strain in arterial wall are small compared with the normal components) and considering the volumetric incompressibility of the vessel, we can obtain the principal stretch ratios λ as (19): (2) where *R*, *Z*, *r*, and *z* are the radial and axial coordinates (the origin is chosen at one end of the vessel) at zero-stress and loaded states, respectively (Fig. 1), and subscript θ denotes circumferential direction. The *r* is related to *R* and *r*_{i} (inner radius of the loaded artery) by (3) The nontrivial Green strain *E* components are (4)

It is commonly assumed that the hydrostatic stress of soft tissues behaves purely elastically, and the deviatoric stresses are viscoelastic (3, 12, 26). Following this axiom approximately (because in general the trace of *S*_{jk} is not zero), we write the quasi-linear viscoelastic model (10) as: (5) where **x** (*r*, θ, *z*) are spatial coordinates, *t* is time, *S*_{jk} is the second Piola-Kirchhoff stress, *g*(*t*) is the reduced relaxation function, *S*_{jk}^{e} denotes the elastic stress (function of current strain and implicitly depends on time), which can usually be obtained from an elastic strain energy function *W*(*E*_{jk}) by (6)

In this study, we employ the widely used Fung-type exponential form of strain energy (8), which, after neglecting the shear components, can be written as: (7) *C* and *b*s are material constants. Without loss of generality, we choose a relatively simple stress relaxation function (8) where τ is the relaxation time, and the coefficients *g*_{0} and *g*_{1} satisfy *g*_{0} + *g*_{1} = 1.

Cauchy stresses σ (when shear deformation is absent) can be computed from the second Piola-Kirchhoff stresses as (8, 19): (9) where *H*(*r*) is an arbitrary scalar that must be determined from boundary conditions.

It should be emphasized that we only model the passive behavior of arteries in this work. Hence, the active stress from the smooth muscle cells (e.g., Ref. 19) is not included in *Eq. 9*. The most general form of a quasi-linear viscoelastic model should be written with a tensorial relaxation function (10). Here, we have simplified it with an isotropic relaxation function (*Eqs. 5* and *8*). In principle, these issues can be taken into account without major difficulties.

Due to the viscoelastic behavior of arteries loaded with pulsatile blood pressure, the radius *r* and axial stretch λ_{z} are functions of time. Considering that the in vivo axial stretch ratio should not change significantly in the cardiac cycle, λ_{z} is assumed to be constant (7, 13). Therefore, the radius *r*(*t*) is determined merely by the inner radius *r*_{i}(*t*), as shown in *Eq. 3*. The radial velocity *v*_{r} and acceleration *a*_{r} of the arterial wall are defined by: (10) The equation of motion in the radial direction can be written as (7, 13): (11) where ρ denotes the mass density of the arterial wall. *Equations 3*–*11* can be solved numerically with pressure boundary conditions on the inner and outer wall surfaces. The detailed numerical procedure is provided in the appendix.

#### Simulation parameters.

The data for rabbit thoracic artery reported by Chuong and Fung (8) were chosen in our analysis: *L*_{i} = 9.75 mm, *L*_{o} = 11.25 mm, Φ = 108.6° (χ = 2.521), λ_{z} = 1.691, *C* = 28.0 kPa (22.4/*g*_{0}), *b*_{1} = 1.0672, *b*_{2} = 0.4775, *b*_{3} = 0.0499, *b*_{4} = 0.0903, *b*_{5} = 0.0585, and *b*_{6} = 0.0042. It is noted that *C* was scaled to make the product *g*_{0}*C* = 22.4 kPa to ensure that the fully relaxed moduli of the viscoelastic model correspond to the elastic moduli obtained in Chuong and Fung for the preconditioned artery. The rationalization for this choice can be found in the appendix.

The mass density ρ was taken to be 0.9 g/cm^{3}, which has been used by Chaudhry et al. (7). Since viscoelastic properties for rabbit thoracic artery in Ref. 8 are not available, typical parameters of the reduced stress relaxation function were selected as τ = 0.1 s, *g*_{0} = 0.8, and *g*_{1} = 0.2 (27, 12). The *g*_{0} = 0.8 is a typical value for the porcine coronary artery reported by Veress et al. (27), the canine aorta reported by Tanaka and Fung (24), and the human coronary arteries considered by Holzapfel et al. (12). Other parameters, such as τ = 0 s [*Eq. 8* reduces to *g*(*t*) = *g*_{0}] and τ = 1.0 s were also employed to study the influence of relaxation time on the hysteresis curves.

To mimic a realistic pulsatile blood pressure, we considered a wave pulse used by Humphrey and Na (13), which is a Fourier series given by: (12) where *A*_{1} = −0.0574, *A*_{2} = −0.0504, *A*_{3} = −0.0105, *A*_{4} = −0.0063, *A*_{5} = −0.0022, *B*_{1} = 0.0668, *B*_{2} = 0.0095, *B*_{3} = −0.0114, *B*_{4} = −0.0032, *B*_{5} = −0.0018, ω = 7.75/s, and *t*_{0} = 0.42 s (the last term is zero in Ref. 13).

The above waveform was used to generate the following pulsatile blood pressure P, which was applied to the inner arterial wall surface as boundary condition: (13) where P_{ib} is the base level of intravascular pressure, and P_{ia} is a scaling factor for the amplitude of the pulse.

All arteries are more or less constrained by perivascular tethering at the in vivo state. Humphrey and Na (13) considered an exponential function to model the tethering, which results in a pulsatile external pressure on the arterial wall. The passive perivascular tethering can be approximated by the same wave pulse as the internal pressure, with a different degree of compression, namely, (14) where P_{ob} and P_{oa} are the base level and the amplitude of the external pressure, respectively.

We chose P_{ib} = 100.6, P_{ia} = 161.2, P_{ob} = 35.1, and P_{oa} = 24.2 (in units of mmHg, 1 mmHg = 0.1333 kPa), which yield an inner pressure of 80–120 mmHg and an outer pressure of 32–38 mmHg (Fig. 2). This external constraint is an average of the tethering considered by Humphrey and Na (13).

## RESULTS

Using the above material parameters and boundary conditions, we first studied the hysteresis behavior of the arterial wall. Figure 3 shows the plots of inner radius and inner pressure (outer pressure also changes, see Fig. 2) for relaxation time τ = 0, 0.1, and 1.0 s. It is seen that the area in the hysteresis loop (an indicator for phase delay between stress and strain) is zero for τ = 0 s and is very small for τ = 1.0 s. In these two cases, the arterial wall behaves as a purely elastic or a nearly elastic material, respectively. In the case of τ = 0.1 s, a larger hysteresis area is observed, implying that the viscoelastic response is present. This phenomenon will be discussed later.

Figure 4*A* shows the temporal evolution of Green strains at the inner surface for τ = 0.1 s, and Fig. 4*B* shows the corresponding pulsatile Cauchy stresses. It is noted that the axial (*z*) stress as well as the circumferential (θ) stress and strain change dynamically with the pulsatile pressure boundary conditions. The dynamic stress and strain are not in phase when viscoelasticity is present, as shown in Fig. 4*C*, where the circumferential stress-strain curves for τ = 0 and 0.1 s are plotted. Note that the magnitudes of stress and strain as well as their variations (difference between the maximum and minimum values) for τ = 0.1 s are smaller than those for τ = 0 s. This implies that the peak stress induced by a viscoelastic response is smaller than that caused by an elastic response at physiological loading. Consequently, the wall stress and strain could be chronically reduced during the course of pulsatile loading, if the viscoelastic property is included.

The energy per unit volume dissipated in one loading cycle *W*_{d} (which is related to heat production rate in vessel wall) can be computed by integrating the second Piola-Kirchhoff stress *S*_{jk} with respect to Green strain *E*_{jk}: (15) where the last loading pulse is chosen (i.e., *t*_{1} = 2.03 s and *t*_{2} = 2.83 s). We obtain *W*_{d} ≈ 0 and *W*_{d} = 0.37 kPa for cases τ = 0 and 0.1 s, respectively. This indicates that finite energy is dissipated into the vessel wall in the viscoelastic model (τ = 0.1 s), but not in the elastic model (τ = 0 s). Figure 5 shows the corresponding plot of circumferential second Piola-Kirchhoff stress and Green strain for τ = 0 and 0.1 s. The shape of these curves is similar to that in Fig. 4*C*, with different stress magnitudes.

To further elucidate the viscoelastic and elastic responses, temporal evolutions of radial velocities and accelerations at the inner surface of the arterial wall for τ = 0 and 0.1 s are plotted in Fig. 6, *A* and *B*, respectively. Consistent with the finding of Humphrey and Na (13), the inertial force in the arterial wall is negligible. The maximum radial acceleration is smaller than 0.92 m/s^{2} (Fig. 6*B*), the dynamic stress caused by inertia is on the order of ρ*a*_{r}*r*_{m} ≈ 2 Pa (*r*_{m} = 2.37 mm is a typical radius of the artery, see Fig. 3), which is insignificant compared with the wall stresses (Fig. 4*B*).

## DISCUSSION

It is well known that morphological changes of large arteries (remodeling) are related to fatigue of load-bearing elastin fibers in cyclic loading (see review in Ref. 7). Based on a comparison between atherosclerotic plaque rupture and engineering fatigue phenomenon, Bank et al. (2) hypothesized that fatigue is a critical mechanism of plaque rupture. Finite-element analysis of an atherosclerotic artery supports the fatigue failure hypothesis (28).

The fatigue of arteries is a chronic failure process caused by numerous cyclic blood pressures. At normal human heart rate, the arterial wall experiences ∼37 millions of pulsatile loading cycles annually. Versluis et al. (28) found that the rate of arterial crack growth depends on both the mean and pulse pressure. They modeled plaque crack growth per loading cycle with the Paris relation. According to the formulation and model parameters in Versluis et al., the crack growth rate ratio (if a crack exists) can be estimated as: (16) where Δ*c*_{v} and Δ*c*_{e} are crack length growth per loading cycle for the viscoelastic (τ = 0.1 s) and elastic (τ = 0 s) responses, respectively; σ_{max}^{v} = 140.3, σ_{min}^{v} = 61.7, σ_{max}^{e} = 141.9, and σ_{min}^{e} = 60.9 (in units of kPa) are circumferential stresses [maximum (max), minimum (min), viscoelastic (*v*), elastic (*e*), respectively] at the inner surface.

It is noted that stresses and strains are largest at the inner surface; e.g., the circumferential stresses (kPa) are σ_{max}^{v} = 109.5, σ_{min}^{v} = 53.9, σ_{max}^{e} = 110.5, and σ_{min}^{e} = 53.3 in the middle of the wall (Δ*c*_{v}/Δ*c*_{e} ≈ 95%), σ_{max}^{v} = 90.8, σ_{min}^{v} = 48.8, σ_{max}^{e} = 91.5, and σ_{min}^{e} = 48.3 at the outer surface (Δ*c*_{v}/Δ*c*_{e} ≈ 96%). In view of the transmural stress gradient, cracks or other damages will likely occur first on the inner side of the arterial wall.

Although crack growth rate for τ = 0.1 s (viscoelastic case) is only 5% less than that for τ = 0 s (elastic case), the arterial fatigue life can be significantly extended, if viscoelasticity is present, since the fatigue failure depends on stress nonlinearly and chronically (2, 28). The accumulative difference in damage due to fatigue can be significant after millions of loading cycles in elastic and viscoelastic arteries. Because stress concentration occurs at the crack tip, the above estimate based on stress distribution without a crack is likely to underestimate the difference. Moreover, the influence of viscoelasticity may be larger than the analyzed case. The reduced relaxation function for subcutaneous tissue of rat decreases to *g*_{0} = 0.6 after 60 s (15), and *g*_{0} = 0.7 is possible after 15 s for porcine coronary artery (27). The dissipated energy in porcine aortic and carotid arteries can be larger than the stored energy (21), which implies that *g*_{0} may have a smaller value than the selected 0.8. A parameter sensitivity analysis shows that Δ*c*_{v}/Δ*c*_{e} ≈ 91% (with σ_{max}^{v} = 139.2 kPa and σ_{min}^{v} = 62.2 kPa) when *g*_{0} = 0.7 and Δ*c*_{v}/Δ*c*_{e} ≈ 88% (with σ_{max}^{v} = 138.0 kPa and σ_{min}^{v} = 62.7 kPa) when *g*_{0} = 0.6. In such cases, the viscoelastic effect is more substantial.

Despite the filtering effect in an instantaneous change of loading thought to prevent mechanical failure (1, 21), the viscoelasticity has not been related to the fatigue of arteries. Based on the present findings, we hypothesize that the reduction of stress and strain due to viscoelastic property of the arterial wall is beneficial for the fatigue life of arteries. This hypothesis can be generalized to other living tissues; i.e., viscoelasticity is present in biological tissues to reduce fatigue failure. Our rationale is that the damage accumulation caused by fatigue could be significantly smaller in viscoelastic response than in purely elastic response after numerous loading cycles. This could possibly occur when an initial injury such as atherosclerotic lesion exists in the arterial wall. The stress concentration in the vicinity of the injured region will undoubtedly lead to accelerated damage/crack growth, especially if the diseased vessel loses its viscoelastic property (e.g., hardened arterial wall due to calcification).

Figure 6 shows that a viscoelastic response (τ = 0.1 s) can reduce the vibration of the arterial wall (smaller radial velocity and acceleration). The smaller wall movement may facilitate mass transport between blood and the artery, since the wall vibration may reduce transient time near the wall. Therefore, it is hypothesized that the recovery of an injured vessel (which needs nutrition from blood) may be improved by viscoelasticity. The implications of the present study to fatigue life of blood vessels, however, are speculative. The validity of the proposed hypothesis must be studied, either by experimental measurements or by an analysis considering cracks, which is beyond the scope of the current model.

One limitation of the proposed model is that the rate-insensitive viscoelastic behavior is not taken into account. In principle, this feature can also be incorporated by either a generalized Maxwell model that includes various relaxation times (11, 12) or by a continuous spectrum relaxation function (3, 10). Since we have invoked the viscoelastic behavior in the relaxation function with a single relaxation time (*Eq. 8*), the salient features (hysteresis behavior and stress reduction) will not change, if a frequency-insensitive viscoelastic model is employed.

Another limitation of the present model is that it does not account for continual turnover, repair, regeneration, and remodeling of the vessel wall components. Thus the results should be interpreted based on the assumption that the cellular and metabolic activities are identical in the viscoelastic and elastic models. The present analysis provides a purely mechanical perspective based on the assumption that the arterial wall (8) has passive viscoelastic properties in *Eq. 5*.

In summary, a quasi-linear viscoelastic model was used to investigate the dynamic stress and strain in an arterial wall undergoing pulsatile blood pressure loading. The numerical results demonstrate that the magnitudes of physiological wall stress and strain can be reduced, if the wall is viscoelastic. This implies that viscoelasticity in biological tissues may be beneficial for the fatigue life. These findings may have important implications for mechanical homeostasis, which undoubtedly affects growth and remodeling when the mechanical environment is perturbed (16).

## APPENDIX

### Numerical Procedure

Substituting the reduced relaxation function (*Eq. 8*) into the viscoelastic model (*Eq. 5*) and assuming that *t* = −∞ corresponds to zero-stress state, we obtain (A1) where the first term on the right-hand side *g*_{0}*S*_{jk}^{e} is an elastic response, and the second term *g*_{1}*T*_{jk} is a viscous response, which vanishes at full relaxation. The viscous term is (A2)

We compute the viscous stress incrementally with a procedure similar to those in Simo (22) and Vena et al. (26). Specifically, we use the stresses at time *t* to compute the stresses at time *t* + Δ*t*. Keeping in mind that the elastic stress *S*_{jk}^{e} (*t* + Δ*t*)can be easily obtained from the strain energy function (*Eq. 6*), we only need to compute the viscous stress. In particular, (A3) The first integral can be readily written as *e*^{−Δt/τ} *T*_{jk}(*t*), and the second can be estimated by (26): (A4) where (A5) The integration in *Eq. A4* can be obtained easily. Finally, the viscous stress in *Eq. A3* becomes (A6) Using *Eq. A6*, the second Piola-Kirchhoff stress *S*_{jk}(*t* + Δ*t*) in *Eq. A1* can be calculated.

Similar to the results in Rachev and Hayashi (19), the current Cauchy stresses (without smooth muscle activity) can be obtained when taking into account the boundary condition σ_{rr}(*r*_{i}) = −P_{i}(*t*) in *Eqs. 11* and *9*: (A7) (A8) (A9) According to *Eqs. A7*–*A9*, the Cauchy stress components at time *t* + Δ*t* are attainable if the corresponding *S*_{jk} and *a*_{r} are known (note that P_{i} is a prescribed boundary condition).

To compute the radial acceleration *a*_{r}, Chaudhry et al. (7) and Humphrey and Na (13) employed the semi-inverse method developed by Demiray and Vito (9), in which the dynamic inner radius of the artery is approximated by a Fourier expansion with given coefficients. Here, we exploit the Newmark-beta method (18) in the numerical integration algorithm for velocity and displacement. When parameter β = 1/4 (unconditionally stable scheme) is selected, the Newmark-beta formulation becomes (A10) (A11)

Since the artery is assumed to be incompressible, the radius *r*(*t*) at an arbitrary material point along the wall thickness can be calculated with *Eq. 3*, if the inner radius *r*_{i}(*t*) is known. Our numerical procedure can be described as follows:

) Discretize arterial wall thickness

*R*_{o}−*R*_{i}into*N*nodes (*N*= 21).) Guess

*r*_{i}(*t*+ Δ*t*) and compute*r*_{n}(*t*+ Δ*t*) (*n*= 1, …,*N*) using*Eq. 3*.) Calculate

*a*_{r}(*t*+ Δ*t*) (*Eq. A11*) and then*v*_{r}(*t*+ Δ*t*) (*Eq. A10*).) Compute λ

_{θ}, λ_{z}, λ_{r}(*Eq. 2*) and Green strains (*Eq. 4*) at*t*+ Δ*t*.) Compute σ

_{rr}, σ_{θθ}, and σ_{zz}at time*t*+ Δ*t*(*Eqs. A7*–*A9*).) If 0 ≤ σ

_{rr}(*r*_{o}) + P_{o}(*t*+ Δ*t*) ≤ 10^{−9}kPa [i.e., boundary condition σ_{rr}(*r*_{o}) = −P_{o}(*t*+ Δ*t*) is approximately satisfied, note that*r*_{o}=*r*_{N}], output result and go back to*step 2*for a new time increment. Otherwise repeat*steps 2–7*for the current time.

It is pointed out that, at *t* = 0, we have assumed *T*_{jk} = 0, *v*_{r} = 0, and *a*_{r} = 0, which is equivalent to a fully relaxed static state (dynamic behavior starts from *t* = 0). The arterial wall stress was computed as (*Eqs. 5* and *A1*): (A12) Hence, the choice for *C* was 28.0 kPa (*g*_{0}*C* = 22.4 kPa). As shown in Fig. 2, *a*_{r}(0) = 0 and *v*_{r}(0) = 0 are reasonable assumptions, as the wave pulse has been appropriately shifted by *t*_{0} in *Eq. 12* to yield nearly zero derivatives of load with respect to time at *t* = 0.

## GRANTS

This research was supported in part by National Heart, Lung, and Blood Institute Grant 2 R01 HL-055554-11.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society