**1.2. Clinical significance of blood vessel stiffness**

Arterial stiffness indexes, such as PWV (pulse wave velocity: velocity of pressure wave propagation in circulatory systems), PP (pulse pressure: systolic blood pressure minus diastolic blood pressure), and AIx (augmentation index: index showing the effect of reflecting waves coming from peripheries) have received a lot of attention because they can indicate the risks for cardiovascular diseases (Oliver & Webb, 2003). These clinical indexes, which are obtained by analyzing blood pressure pulse waveforms *in vivo,* are used for predicting the prognosis of cardiovascular diseases and thus analyses of pulse waveform are clinically important. In order to understand the basis of these indexes, it is necessary to understand how changes in mechanical properties, such as blood vessel viscoelasticity, can affect pulse waveforms in the systemic arteries.

### **1.3. Research purpose and contents of this chapter**

We thereby proposed a nonlinear one-dimensional numerical simulation model which can accurately calculate the viscous resistance of unsteady flow and the viscoelasticity of the tube wall. We have shown that the numerical model can describe well wave propagation in silicone tubes representing blood vessels (Kitawaki et al., 2003). The result showed that the viscoelasticity of the vessel wall plays an important role in the form of a pulsatile wave (Kitawaki & Shimizu, 2005); therefore, it is necessity to consider the viscoelastic effect accurately in the quantitative investigation of changes in pulsatile waves *in vivo*.

In **section 2**, we construct a one-dimensional numerical simulation model that takes into account unsteady viscosity and the generalized viscoelastic model from the theory of the fluid equation (Kitawaki et al., 2003; Kitawaki & Shimizu, 2006). In addition, the high-speed calculation method is described. When the wall of a blood vessel deforms finitely due to changes in the internal pressure, the wall's physical properties, such as deformation compliance and viscoelasticity, change nonlinearly (Hayashi et al., 1980; Sato & Ohshima, 1985). In **section 3**, it was checked whether our one-dimensional model would be able to simulate the pulse wave propagation of small pressure waves in silicone tubes even when their deformation compliance and viscoelasticity changed independently, using appropriate values of the viscoelastic parameter of the silicone rubber tube for numerical simulation (Kitawaki & Shimizu, 2006). In **section 4**, the influence of viscoelasticity change on periodic pulsatile wave propagation was studied (Kitawaki & Shimizu, 2009)**.** The purpose of the section was to investigate the effect of viscoelasticity change on periodic pulsatile wave quantitatively. For this purpose, we studied the pulse wave propagation of periodic small pressure waves using a silicone tube connected with terminal resistance, and obtained the waveform changes of pulsatile waves due to the change of mechanical properties, including the viscoelasticity of the tube. We then compared the experimental results with numerically calculated results using a one-dimensional numerical simulation model with terminal resistance treatment. Finally, we studied the effect of vessel wall viscoelasticity on the propagation of a periodic pulsatile wave by comparing numerical simulation results between the difference of viscoelastic models and viscoelastic parameters. The final **section 5**, describes the conclusion of this chapter.

#### **2. Theory**

188 Viscoelasticity – From Theory to Biological Applications

**1.2. Clinical significance of blood vessel stiffness** 

affect pulse waveforms in the systemic arteries.

**1.3. Research purpose and contents of this chapter** 

Consequently, in order to construct a numerical simulation model of intravascular flow for quantitative analysis with viscoelasticity of the blood vessel wall, it is necessary to use a

Arterial stiffness indexes, such as PWV (pulse wave velocity: velocity of pressure wave propagation in circulatory systems), PP (pulse pressure: systolic blood pressure minus diastolic blood pressure), and AIx (augmentation index: index showing the effect of reflecting waves coming from peripheries) have received a lot of attention because they can indicate the risks for cardiovascular diseases (Oliver & Webb, 2003). These clinical indexes, which are obtained by analyzing blood pressure pulse waveforms *in vivo,* are used for predicting the prognosis of cardiovascular diseases and thus analyses of pulse waveform are clinically important. In order to understand the basis of these indexes, it is necessary to understand how changes in mechanical properties, such as blood vessel viscoelasticity, can

We thereby proposed a nonlinear one-dimensional numerical simulation model which can accurately calculate the viscous resistance of unsteady flow and the viscoelasticity of the tube wall. We have shown that the numerical model can describe well wave propagation in silicone tubes representing blood vessels (Kitawaki et al., 2003). The result showed that the viscoelasticity of the vessel wall plays an important role in the form of a pulsatile wave (Kitawaki & Shimizu, 2005); therefore, it is necessity to consider the viscoelastic effect

In **section 2**, we construct a one-dimensional numerical simulation model that takes into account unsteady viscosity and the generalized viscoelastic model from the theory of the fluid equation (Kitawaki et al., 2003; Kitawaki & Shimizu, 2006). In addition, the high-speed calculation method is described. When the wall of a blood vessel deforms finitely due to changes in the internal pressure, the wall's physical properties, such as deformation compliance and viscoelasticity, change nonlinearly (Hayashi et al., 1980; Sato & Ohshima, 1985). In **section 3**, it was checked whether our one-dimensional model would be able to simulate the pulse wave propagation of small pressure waves in silicone tubes even when their deformation compliance and viscoelasticity changed independently, using appropriate values of the viscoelastic parameter of the silicone rubber tube for numerical simulation (Kitawaki & Shimizu, 2006). In **section 4**, the influence of viscoelasticity change on periodic pulsatile wave propagation was studied (Kitawaki & Shimizu, 2009)**.** The purpose of the section was to investigate the effect of viscoelasticity change on periodic pulsatile wave quantitatively. For this purpose, we studied the pulse wave propagation of periodic small pressure waves using a silicone tube connected with terminal resistance, and obtained the waveform changes of pulsatile waves due to the change of mechanical properties, including the viscoelasticity of the tube. We then compared the experimental results with numerically

accurately in the quantitative investigation of changes in pulsatile waves *in vivo*.

nonlinear distributed parameter model to be able to include the viscoelasticity.

#### **2.1. Basic equations**

When constructing the numerical model, we neglect the effects of bends of vessels. We also assume that the tube does not leak and that the flow is axisymmetric and incompressible. Under these conditions, the equations of continuity and momentum conservation of the onedimensional model are given by (Olufsen, 1999),

$$\frac{\partial A}{\partial t} + \frac{\partial Q}{\partial \mathbf{x}} = 0 \tag{1}$$

$$
\frac{\partial \mathbb{Q}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left( \frac{\mathbb{Q}^2}{A} \right) + \frac{A}{\rho} \frac{\partial p}{\partial \mathbf{x}} + F\_t = \mathbf{0} \tag{2}
$$

where *A* is the cross-section of the tube, *Q* is the mean sectional flow volume, *p* is the mean pressure, *t* is time, *x* is distance along the vessel axis, *ρ* is the fluid density, and *Ft* is the viscous resistance.

By assuming that we are dealing with a Newtonian fluid, and an oscillating flow velocity distribution in a cylindrical tube, using the Womersley model, the viscous resistance *Ft* is given by (Zielke, 1968),

$$F\_t = 4\pi\nu \left\{ 2V + \int\_0^t \frac{\partial V}{\partial t} W(t - \tau) d\tau \right\} \tag{3}$$

where *V* is the mean sectional velocity, *ν* is the kinematic viscosity, and *W(t)* is a weight function. In a rigid cylindrical tube,

$$W(t) = \lim\_{\alpha \to \infty} \int\_{-\pi}^{\pi} \frac{1}{F\_f(\alpha i^{\sqrt{2}}) - 2} e^{i\alpha t} d\alpha \tag{4}$$

where 0 1 () () () *<sup>J</sup> F z zJ z J z* , and *J0*, *J1* are 0th and 1st order Bessel functions of the complex number *z*. *ω*(=2π*f*) is the angular frequency of the flow oscillating at frequency *f*, *R* is the Womersley number and *R* is the tube radius.

For a long wavelength, the flow velocity distribution in a distensible tube is similar to that in a rigid tube. Therefore, the weight function in a distensible tube can be approximated by Eq. (4),

regardless of changes in the tube's cross-section induced by the internal pressure. Note that the radius expressed as the Womersley number becomes a function of position and time, and consequently the weight function itself also becomes a function of position and time.

#### **2.2. Deformation models of a wall**

The tube law that describes the relationship between the tube cross-section and the internal pressure can be described by an elastic model and two viscoelastic models as below. In actuality, for a silicone tube, the relationship varies with internal pressure change, and shows nonlinearity. However, by assuming small local transformations, the following models can be used to approximate a linear model.

#### *2.2.1. Elastic model*

When the tube is perfectly elastic, the tube law of the elastic model can be expressed by,

$$
\Delta p = \frac{1}{Cs} \Delta A \tag{5}
$$

Numerical Simulation Model with Viscoelasticity of Arterial Wall 191

In order to obtain the tube law in one-dimensional flow analysis, the tube law of the generalized viscoelastic model is derived from the complex viscoelastic coefficient as follows. The complex viscoelastic coefficient of the generalized viscoelastic model *E(s)* can

> <sup>0</sup> 1 1 ( ) (1 ' ) (1 ) *n n i i i i*

where *E0* is the static elastic coefficient, *s* is angular frequency (=*iω*), and *τi* and, *τ'i* are relaxation time parameters which express the viscoelasticity of the tube wall. This equation

> 1 () 1 1 / *n i i i*

0 1 () () 1 <sup>2</sup> 1 /

where *h* is the tube wall thickness, *Δp(s)* and *ΔA(s)* are Laplace transformations of the pressure and cross-section from a reference value. The compliance of the tube deformation

By transforming Eq. (10) back to the time domain by inverse Laplace transformation, we

<sup>0</sup> <sup>1</sup> 1 ( ) *<sup>i</sup> <sup>n</sup> <sup>t</sup> t u i i dA u p A fe du Cs dt*

The following high-speed calculation method that was obtained for the viscous resistance term of rigid tube (Kagawa et al., 1983) was applied to the convolution integrals which appeared in the viscous resistance term in Eq. (3) and the viscoelastic terms in Eq. (7). By approximating the

( )/

 

*hE <sup>s</sup> ps As f RA <sup>s</sup>*

*s*

 

> *n i i i*

 (10)

(11)

0 ( ) *<sup>i</sup> <sup>k</sup> <sup>n</sup> i i W t me*

; the

(8)

(9)

*Es E s s* 

0

.

0

*<sup>s</sup> Es E f*

**2.3. Generalized viscoelastic model** 

can be transformed:

where

*f*

is <sup>0</sup>

2 *s RA <sup>C</sup> hE* .

0

be expressed by (Westerhof and Noordergraaf, 1970),

1 1 ( ') ( )

*i k*

Using this formula, the tube law in the frequency domain becomes

obtain the tube law of the generalized viscoelastic model.

weight function in Eq. (4) by the sum of the exponential function,

viscous resistance in Eq. (3) can be expressed by recursive formulations as follows,

**2.4. High-speed calculation method** 

 

*n n i i k i ik k k*

where Δ*p* and Δ*A* are the differences in the pressure and cross-section relative to reference values, and 0 *<sup>A</sup> Cs p A* is tube deformation compliance at reference cross-section *A0*.

#### *2.2.2. Voigt model*

By assuming that the viscoelasticity of the tube causes a phase lag between the applied pressure and resulting change in cross-section of the tube, from the Voigt model, the tube deformation law can be expressed by,

$$
\Delta p = \frac{1}{\text{Cs}} \left\{ \Delta A + \tau\_V \frac{\partial A}{\partial t} \right\} \tag{6}
$$

where *τV* is the relaxation time that accounts for the phase lag.

#### *2.2.3. Generalized viscoelastic model*

For complex viscoelasticity, as is the case for blood vessels, the generalized viscoelastic model can be applied (Westerhof and Noordergraaf, 1970). The following tube law of the generalized viscoelastic model can be derived, as shown in the next paragraph.

$$\Delta p = \frac{1}{\text{Cs}} \left| \Delta A + \int\_0^t \sum\_{i=1}^n f\_i e^{-(t-u)/\tau\_i} \frac{dA(u)}{dt} du \right| \tag{7}$$

The viscoelastic property of the tube wall is reflected in the second term of Eq. (7), which contains both the dynamic viscoelasticity parameter *fi* and the relaxation time parameter *τi*.

#### **2.3. Generalized viscoelastic model**

190 Viscoelasticity – From Theory to Biological Applications

**2.2. Deformation models of a wall** 

*2.2.1. Elastic model* 

values, and

*2.2.2. Voigt model* 

models can be used to approximate a linear model.

0

*p A*

deformation law can be expressed by,

*2.2.3. Generalized viscoelastic model* 

*<sup>A</sup> Cs*

regardless of changes in the tube's cross-section induced by the internal pressure. Note that the radius expressed as the Womersley number becomes a function of position and time, and

The tube law that describes the relationship between the tube cross-section and the internal pressure can be described by an elastic model and two viscoelastic models as below. In actuality, for a silicone tube, the relationship varies with internal pressure change, and shows nonlinearity. However, by assuming small local transformations, the following

When the tube is perfectly elastic, the tube law of the elastic model can be expressed by,

<sup>1</sup> *p A Cs*

where Δ*p* and Δ*A* are the differences in the pressure and cross-section relative to reference

By assuming that the viscoelasticity of the tube causes a phase lag between the applied pressure and resulting change in cross-section of the tube, from the Voigt model, the tube

> *<sup>A</sup> p A Cs t*

For complex viscoelasticity, as is the case for blood vessels, the generalized viscoelastic model can be applied (Westerhof and Noordergraaf, 1970). The following tube law of the

> <sup>0</sup> <sup>1</sup> 1 ( ) *<sup>i</sup> <sup>n</sup> <sup>t</sup> t u i i dA u p A fe du Cs dt*

The viscoelastic property of the tube wall is reflected in the second term of Eq. (7), which contains both the dynamic viscoelasticity parameter *fi* and the relaxation time parameter *τi*.

1

generalized viscoelastic model can be derived, as shown in the next paragraph.

where *τV* is the relaxation time that accounts for the phase lag.

is tube deformation compliance at reference cross-section *A0*.

*V*

( )/

(7)

 

(5)

(6)

consequently the weight function itself also becomes a function of position and time.

In order to obtain the tube law in one-dimensional flow analysis, the tube law of the generalized viscoelastic model is derived from the complex viscoelastic coefficient as follows. The complex viscoelastic coefficient of the generalized viscoelastic model *E(s)* can be expressed by (Westerhof and Noordergraaf, 1970),

$$E(\mathbf{s}) = E\_0 \prod\_{i=1}^{n} (1 + s\tau\_i') \bigg/ \prod\_{i=1}^{n} (1 + s\tau\_i) \tag{8}$$

where *E0* is the static elastic coefficient, *s* is angular frequency (=*iω*), and *τi* and, *τ'i* are relaxation time parameters which express the viscoelasticity of the tube wall. This equation can be transformed:

$$E(\mathbf{s}) = E\_0 \left( 1 + \sum\_{i=1}^{n} f\_i \frac{\mathbf{s}}{\mathbf{s} + \mathbf{1} / \tau\_i} \right) \tag{9}$$

where 1 1 ( ') ( ) *n n i i k i ik k k i k f* .

Using this formula, the tube law in the frequency domain becomes

$$
\Delta p(\text{s}) = \frac{hE\_0}{2RA\_0} \cdot \Delta A(\text{s}) \cdot \left\{ 1 + \sum\_{i=1}^n f\_i \frac{\text{s}}{\text{s} + 1/\text{ }\tau\_i} \right\} \tag{10}
$$

where *h* is the tube wall thickness, *Δp(s)* and *ΔA(s)* are Laplace transformations of the pressure and cross-section from a reference value. The compliance of the tube deformation is <sup>0</sup> 2 *s RA <sup>C</sup> hE* .

By transforming Eq. (10) back to the time domain by inverse Laplace transformation, we obtain the tube law of the generalized viscoelastic model.

$$
\Delta p = \frac{1}{\text{Cs}} \left| \Delta A + \int\_0^t \sum\_{i=1}^n f\_i e^{-(t-u)/\tau\_i} \frac{dA(u)}{dt} du \right| \tag{11}
$$

#### **2.4. High-speed calculation method**

0

The following high-speed calculation method that was obtained for the viscous resistance term of rigid tube (Kagawa et al., 1983) was applied to the convolution integrals which appeared in the viscous resistance term in Eq. (3) and the viscoelastic terms in Eq. (7). By approximating the weight function in Eq. (4) by the sum of the exponential function, 0 ( ) *<sup>i</sup> <sup>k</sup> <sup>n</sup> i i W t me* ; the viscous resistance in Eq. (3) can be expressed by recursive formulations as follows,

$$F\_t(t) = 4\pi\nu \left\{ 2V(t) + \sum\_{i=0}^{k} y\_i(t) \right\} \tag{12}$$

Numerical Simulation Model with Viscoelasticity of Arterial Wall 193

valve

0 10 20 30 40 50 Frequency (Hz)

(c) Frequency component of the movement

Water tank

Figure 2(a) is a schema of the piston pump. The piston pump was driven by a computercontrolled stepping motor, and was capable of generating various waveforms with various flow volumes. The piston cylinder receives the backpressure from the water inside the tube. Therefore, displacement of the inner cylinder was measured by a laser displacement sensor

> Pressure sensor (P1)


Normalized Fourier component

Pressure sensors

(with pressure sensor at the center of the tube)

**P2 P3**

1.45m silicone tube

The tube, piston pump, and water tank were filled with water. Baseline internal pressures were set by adjusting the water head of the tank. The piston pump generated a single impulse. The flow volume at the inlet of the experiment tube was determined as shown in Fig. 2(b); when the impulse time *ti* was set to 0.1 s, and the maximum flow volume *Qm* was 5.0 ml/s. As seen in Fig. 2(c) which shows the Fourier transformation of the flow volume when a single impulse was generated, the highest frequency component of the flow volume was 30 Hz. Signals from the three pressure sensors and the displacement sensor were recorded for 10 seconds by a PC at a sampling rate of 1 kHz. The trials, that is, generation of a single impulse by the pump, were performed at increasingly higher baseline pressures of 2.5, 5.0, 7.5 and 10.0 kPa at 90-minute intervals. The trials were generated more than 70 minutes after changing the baseline pressures because 30 minutes were necessary for the viscoelasticity to return to a steady state. The baseline pressures were established by

opening the valve, and the valve was closed just before the start of each trial.

**Figure 1.** Experimental apparatus

Piston pump

**P1**

0.05m brass rigid tube

(Keyence, LK-030).

**Figure 2.** Piston pump

0.025 0.075 0.050 0.1

Flow (ml/s)

Stepping Motor (Controlled by PC)

Displacement sensor

> 5.0 2.5

> > Time (s)

(b) movement of piston pump

(a) a cross section of piston pump

2

Inner cylinder

*t t t t Q t t t*

 

2 2 5{8( / 0.1) } ( 0.025) 5{1 2 (2 / 0.1 1) } (0.025 0.075) ( ) 5{8( / 0.1 1) } (0.025 0.1) 0 (0.1 )

barrel

*t*

*Experimental conditions* 

$$\begin{cases} y\_i(t) = 0 & \text{( $t = 0$ )}\\ y\_i(t + \Delta t) = e^{-n\_i \Lambda \tau} y\_i(t) & \text{( $t > 0$ )}\\ \quad + m\_i e^{-n\_i(\Lambda \tau/2)} \{V(t + \Delta t) - V(t)\} \end{cases} \tag{13}$$

where *τ*=*νt/R2* and *Δτ=νΔt/R2* are normalized time and time step, respectively.

We need to determine the value of term number *k* in the Eq. (12) from the value of *Δτ* with consideration of approximation accuracy. In present experimental condition, the term number *k* was 10, because *Δτ* was calculated to be 1.2x10-5.

The weight function of the convolution integral in the tube law of Eq. (7) can be expressed as a summation of the exponential function. Therefore, the tube law of the recursive formulations is as follows,

$$p(t) - p\_0 = \frac{1}{\text{Cs}} \left\{ (A - A\_0) + \sum\_{i=1}^{n} z\_i \right\} \tag{14}$$

$$\begin{cases} z\_i(t) = 0 & \text{( $t = 0$ )}\\ z\_i(t + \Delta t) = e^{-\Lambda t/\tau\_i} z\_i(t) & \text{( $t > 0$ )}\\ \quad + f\_i e^{-\Lambda t/2\tau\_i} \{A(t + \Delta t) - A(t)\} \end{cases} \tag{15}$$

It is possible to determine the term number *n* in Eq. (14) from the viscoelastic characteristics.
