## Abstract

The pulsatility of coronary circulation can be accurately simulated on the basis of the measured branching pattern, vascular geometry, and material properties of the coronary vasculature. A Womersley-type mathematical model is developed to analyze pulsatile blood flow in diastole in the absence of vessel tone in the entire coronary arterial tree on the basis of previously measured morphometric data. The model incorporates a constitutive equation of pressure and cross-section area relation based on our previous experimental data. The formulation enables the prediction of the impedance, the pressure distribution, and the pulsatile flow distribution throughout the entire coronary arterial tree. The model is validated by experimental measurements in six diastolic arrested, vasodilated porcine hearts. The agreement between theory and experiment is excellent. Furthermore, the present pulse wave results at low frequency agree very well with previously published steady-state model. Finally, the phase angle of flow is seen to decrease along the trunk of the major coronary artery and primary branches toward the capillary vessels. This study represents the first, most extensive validated analysis of Womersley-type pulse wave transmission in the entire coronary arterial tree down to the first segment of capillaries. The present model will serve to quantitatively test various hypotheses in the coronary circulation under pulsatile flow conditions.

- pulse wave transmission
- Womersley's method
- impedance
- admittance
- Fourier transform

the blood flow in the coronary arteries is very pulsatile with zero or even reversing (negative) flow in systole. The two major determinants of coronary flow pulsatility are *1*) the pulsatile aortic pressure, and *2*) the myocardial/vascular interaction in systole. To understand coronary blood flow, we must understand each of these two effects. Because the former is simpler, it is the logical starting point. Hence, the objective of the present study is to develop a validated mathematical model of pulsatile blood flow in the entire coronary arterial tree in diastole in the absence of vessel tone. Kassab et al. (17, 18, 20) have previously measured the morphometric data of the coronary vasculature. Recently, a computer model of the entire coronary arterial tree has been developed by Mittal et al. (22) based on measured morphometric data of Kassab et al. (20). This model will serve as the platform for the present systematic pulse pressure and flow wave transmission analysis.

The theory of pulsatile blood flow was initiated by Leonhard Euler in 1775 and Thomas Young in 1808. Later on, numerous attempts have been made at the study of pulsatile blood flow in the vascular system, e.g., Zamir's Womersley-type analysis (5, 6, 33, 34), and Hughes' one-dimensional theory (15). Basically, there are three major approaches to simulate the pulsatile blood flow and pressure waves in the vascular system: *1*) the Womersley's solution of pulsatile blood flow (1, 6, 9–11, 14, 31, 32), *2*) the one-dimensional wave propagation solution (7, 15, 24, 25, 27, 29), and *3*) the three-dimensional solution of pulsatile blood flow (7). The latter two methods require appropriate outflow boundary conditions to simulate the pulsatile blood flow wave. Different outflow boundary conditions have been developed to satisfy the one- or three-dimensional equations, such as the pure resistive load boundary conditions (28), the windkessel model (16, 21), the nonreflecting outlet boundary conditions (7), or the impedance boundary conditions (24, 25, 29). All of the boundary conditions may be classified as lumped-parameter model boundary conditions. Despite the usefulness of lumped models, however, they are generally limited to the global aspects of coronary blood flow and fail to reveal the flow and pressure wave distribution through the vascular system. For example, lumped models cannot be used to predict the significant spatial distribution of coronary blood flow. Womersley's solution (1, 6, 9–11, 14, 31, 32) avoids those shortcomings because it assumes that pressure and flow are decoupled. Hence, this approximation can be applied throughout a complex vascular network such as the entire coronary arterial tree.

Here, we shall use the platform of the entire coronary arterial tree (22) to carry out a Womersley-type numerical analysis of pulse wave transmission. The constitutive equation that relates pressure to cross-sectional area for the several largest orders is based on experimental data (19). The Womersley model is modified to account for the measured constitutive relation of coronary arteries. Previous measurements of coronary wall thickness (13) are also incorporated into the numerical model. The predictions of the mathematical model are compared with the experimental measurements on diastolically arrested, vasodilated porcine hearts. The agreement between theory and experiment is excellent. The pulse wave transmission is analyzed along spatial pathways stemming from the trunk through various branches to the first vessel segments of capillaries. The limitations and implications of theory and experiment are contemplated.

## Glossary

*a,b*- Arbitrary constants of integration in the equations for flow and pressure
*A*- Area of a vessel (cm
^{2}) *A*(*n*)- Initial area of a vessel (cm
^{2}) *c,c*_{0}- Wave velocity with/without the viscosity effect
*c*= ·*c*_{0}(*c*_{0}= ) (cm/s) *D*- Diameter of a vessel (cm)
*D*(*n*)- Initial diameter of a given vessel (cm)
*E*(*E*_{stat})- Young's modulus (static) (dyn/cm
^{2}) *E*_{dyn}- Young's modulus (dynamic) (dyn/cm
^{2}) *f*- Wave frequency (Hz)
*F*_{10}*F*_{10}(α) = [2*J*_{1}(*i*^{3/2}α)]/[*i*^{3/2}α*J*_{0}(*i*^{3/2}α)], where*J*_{0}and*J*_{1}are Bessel function of the first kind, and order zero and one, respectively, and α =*R*(*n*) is the Womersley number*h*[*h*(*n*)]- Wall thickness of a vessel
*n*(cm) *L*(*n*)- Initial length of a vessel (cm)
*n*- Vessel segment number
- p
_{0} - External pressure, assumed zero in diastole (mmHg)
- P(
*x*,ω) - Pressure along the axial length
*x*of the vessel (mmHg) - Q(
*x*,ω) - Volumetric flow rate along the axial length of the vessel (ml/min)
*R*- Radius of a vessel (cm)
*R*(*n*)*R*(*n*) =*D*(*n*)/2; Initial radius of a vessel (cm)*t*- Time (s)
*u*_{r}- Radial velocity (cm/s)
*u*_{x}- Axial velocity (cm/s)
*x*- Axial position along the length of a vessel (cm)
*Y*_{0}*Y*_{0}=*A*(*n*)/ρ*c*_{0}; Characteristic admittance (cm^{4}·s/g)*Y*(*x*,ω)*Y*(*x*,ω) = 1/[*Z*(*x*,ω)]; Admittance along the length of a vessel (cm^{4}·s/g)*Z*_{0}- Z
_{0}= 1*Y*0 Characteristic impedance (g·s^{−1}·cm^{−4}) *Z*(*x*,ω)*Z*(*x*,ω) = [P(*x*,ω)]/[Q(*x*,ω)]; Impedance along the length of a vessel (g·s^{−1}·cm^{−4})- α
- α =
*R*(*n*); Womersley number - ρ
- Blood density (g/cm
^{3}) - τ
_{θ} - Mean circumferential stress in a vessel wall (dyn/cm
^{2}) - ω
- Angular frequency (rad/s)
- μ
- Blood viscosity (cp)

## MATHEMATICAL MODEL

#### Governing equations.

The details of mathematical derivations are shown in the appendix. Briefly, the governing equations for flow and pressure may be written as (1) (2) where *a* and *b* are arbitrary constants of integration, *c* = √1 − *F*_{10}(α)·*c*_{0} (*c*_{0} = ) is the wave velocity, *Y*_{0} = *A*(*n*)/ρ*c*_{0} is the characteristic admittance, *Z*_{0} = 1/*Y*_{0} is the characteristic impedance, *Y*_{1} = *Y*_{0}, and *Z*_{1} = *Z*_{0}. Following the transmission line method (TLM) (3), we define the impedance and admittance as: (3) (4) In a given vessel segment, at *x* = 0 and *x* = *L*, we have the following inlet and outlet impedance: (5) and (6) If we combine *Eqs. 5* and *6*, we obtain: (7) *Equation 7* was used to calculate the impedance/admittance in the entire coronary tree from inlet to the capillary vessels.

#### Method of solution.

The characteristic impedance, characteristic admittance, and velocity (including the viscous effect) were first calculated for every vessel segment in the entire coronary arterial tree. There are two or more vessels that emanate from the *j*th junction point (e.g., junction *point B* in Fig. 1) anywhere in the current tree structure. We assume that mass is conserved at each junction and pressure is continuous at the junction, which may be written as: (8) (9) From *Eqs. 8* and *9*, we may write: (10) Once the terminal impedance/admittance of the first capillary is computed, we proceed backward to iteratively calculate the impedance/admittance in the entire coronary tree by using *Eq. 7* (obtained from *Eq. 3*) and *Eq. 10* as demonstrated in Fig. 1*B*. The sawtooth pulsatile pressure produced experimentally by the piston pump was discretized by a Fourier transformation to determine the constants *a* and *b* in *Eqs. 1* and *2*. The flow and pressure were then calculated by using *Eqs. 1* and *2*.

The morphometric and mechanical parameters used for the computation of the pulsatile flow through the arterial tree were adopted from our previous publications (13, 19, 22). The length and diameter of every vessel segment were assigned to the entire coronary arterial tree (22), on the basis of detailed morphometric data. The blood flow density (*ρ*) in the coronary arterial tree was assumed to be 1.23 g/cm^{3} (8). The variation of viscosity with vessel diameter and hematocrit was based on the viscosity model of Pries et al. (26). The coronary wall thickness for every order was adopted from Guo and Kassab's data (13). The static Young's modulus was calculated as ∼7.0 × 10^{6} (dyn/cm^{2}) on the basis of previous experimental data (19) (see details in the appendix). The dynamic Young's modulus was also considered (see details in the appendix). The outlet impedances at the first capillary vessel segments (*order 0*) were computed by the steady value, i.e., 128 μ_{capillary}*L*_{capillary}/π*D*(g·s^{−1}·cm^{−4}). Unless otherwise stated, all computations used previously measured physical properties and parameters as described above.

## EXPERIMENTAL METHODS

#### Animal preparation.

Studies were performed on six 3- to 4-mo-old farm pigs weighing 28–38 kg. The experimental procedures of the animal preparation were similar to those described by Kassab et al. (20). All animal experiments were performed in accordance with national and local ethical guidelines, including the Institute of Laboratory Animal Research guide, Public Health Service policy, Animal Welfare Act, and an approved University of California, Irvine-Institutional Animal Care and Use Protocol.

#### Measurement of pulsatile flow.

After the heart was isolated, it was placed in a cold (0°C) saline bath as shown in Fig. 2Fig. 3. Therefore, we used the piston pump to simulate the in vivo aortic pressure. In the present study, the wave amplitude of the pulsatile pressure was in the range of 50–100 mmHg with a mean pressure of ∼75 mmHg.

To mimic the profile and range of the in vivo pressure wave, we selected tubings of appropriate resistance and compliance arranged between the piston pump and pressure transducer. It was found that the phase difference is very small (at most 20 ms), which was within the range of relative error of the flow transducers.

## RESULTS

The results are divided into three subsections. First, low-frequency flow (ω → 0) is compared with steady-state flow. Second, the pressure produced by the piston pump is applied as the inlet pressure boundary condition to calculate the pressure and flow distribution in the entire coronary arterial tree by using the current mathematical model in conjunction with a Fourier transform. The mathematical predictions are in turn compared with the experimental results. Finally, the effects of various parameters (e.g., frequency, viscosity, etc.) on pulsatile flow are examined for the entire coronary arterial tree.

#### Low-frequency flow (ω → 0) compared with previous steady-state flow model.

It is desirable to verify the computer code for the values of the pressure and flow against available analytical or numerical solutions. In the present study, the low-frequency (ω → 0) pressure and flow were calculated by using the current mathematical model and compared with a previous steady-state model developed by our group (23). The two models were implemented with the same mean inlet pressure of 100 mmHg. The outlet pressures in the steady-state computation were set at 16 mmHg and the outlet impedances in the pulsatile computation were evaluated as 128 μ_{capillary}*L*_{capillary}/π*D*(g·s^{−1}·cm^{−4}) at the first capillary vessel segments (*order 0*). The wave frequency (*f*) of the pulsatile blood flow was assumed to be 0.001 Hz. The entire asymmetric LAD, LCX, and RCA trees were considered in the comparison of the two models. The direct relationship between mean segment flow |Q|_{mean} (average of flow rates at every order) and mean segment pressure |P|_{mean} {average of 0.5 × [|P(0,ω)| + |P(*L*,ω)|] at every order} in the coronary arteries is depicted in Fig. 4. The present pulsatile flow results at low frequency agree very well with previous steady-state model.

#### Experimental validation of mathematical model.

The piston pump was used to impose the pressure and flow wave at the inlet of the left common coronary artery (LCCA) and RCA. As described in *Measurement of pulsatile flow*, it was established that the profile of the sawtooth pressure produced by the piston pump was similar to that measured in vivo when the piston pump was operated at a frequency of 90 beats/min. The Fourier transform was applied to discretize the real measured pressure wave into the complex periodic Fourier series. The mathematical methods, as described in mathematical model, were then implemented to calculate the flow rate at the inlet of the LCCA and RCA trees based on the detailed morphometric data (20). The viscosity (μ) of the perfusion solution was selected as 1.1 cp to mimic our cardioplegic solution containing albumin. Figure 5, *A* and *B*, shows measured pressures (means ± SD; average and SD of the pressure in 6 hearts) and flow rates (means ± SD) in the waveform at the inlet of LCCA and RCA of the respective hearts. We normalized the flows to the heart weight of 150 g for which the morphometric data were obtained. Figure 6, *A* and *B*, depicts the comparison of flow rates between theoretical prediction and experimental measurement (means ± SD; average and SD of the flow rates in 6 hearts) at the inlet of LCCA and RCA. It is found that the theoretical predictions are within 1 SD of the experimental data.

Figure 7 shows a comparison of numerical results predicted by the present model with those of other models (see Refs. 24 and 25 and Refs. 6, 33, and 34). The formulations of various models are extended to the entire coronary arterial tree and with the same inlet pressure wave as shown in Fig. 5. It is noted that there is a small phase angle difference between the current model and Olufsen's model (24, 25) and a large amplitude difference between the current model and Zamir's formulation (6, 33, 34). Zamir's formulations predict values ∼5 times larger than the current model if it does not consider the viscosity effect (Fig. 7*B*). If the viscosity effect is considered in Zamir's formulation, the disparity with our model is even more pronounced.

#### Effect of various parameters on pulsatile blood flow.

To study the effect of wave frequency, the pressure and flow waves at the inlet of an artery were produced by the piston pump at a heart rate of 90, 120, and 180 beats/min. The flow waveform was calculated by using the measured input pressure and the mathematical model and compared with the measured flow waveform. Figure 8, *A–C*, shows the pressure and flow waveform at the inlet of a LAD artery at different heart rates (90, 120, and 180 beats/min). It is found that the flow waveform reveals similar behavior for experimental and numerical results. To evaluate the effect of pulsatile flow on the mean flow rate with blood viscosity [based on the viscosity model of Pries et al. (26)], Table 1 lists mean pressure (average over 2 min) and mean flow rate (average over 2 min) at the inlet of LAD artery at heart rate of 90, 120, 180 beats/min and steady-state flow rate when the inlet pressure is the same as the mean pressure of pulsatile flow.

To study the effect of the harmonic wave frequency (0 to 3 Hz), changes of the complex flow rate [Q_{inlet}(0,ω)] and impedance [*Z*_{inlet}(0,ω)] at the inlet of LAD artery are shown in Table 2. The complex inlet pressure P_{inlet}(0,ω) was set at 100 mmHg with zero phase angle. Figure 9 shows the change of the amplitude and phase angle of the impedance [*Z*_{inlet}(0,ω)] at the inlet of LAD artery with angular frequency (ω = 2π*f*).

The major advantage of Womersley's method is that it can predict the significant spatial distribution of coronary blood flow. Therefore, the blood flow and pressure in vessel segments along the trunk and the primary branches (branches that arise directly from the trunk) at different harmonic frequencies were investigated for the LAD arterial tree (Figs. 10 and 11). The inlet pressure P_{inlet}(0,ω) was set at 100 mmHg with no phase angle at the inlet of the LAD arterial tree. Figure 10*A* shows a schematic of the trunk and the primary branches, several of which were identified alphabetically (i.e., *branches A–C*). Figure 10*B* shows the relationship between the flow rate |Q(*x,*ω)| and the cumulative length [∑*L*(*n*)] of vessel segments from the root (*order 11*) to the first capillary vessel segment (*order 0*) along the trunk (root-AP-capillary) with harmonic wave frequency of 0.001 Hz and 1.5 Hz. Figure 10*C* shows the corresponding data for the phase angle. Figure 11*A* shows the relationship between the flow rate |Q(*x,*ω)| and the cumulative length of vessel segments through the primary branches (*branches A*, *B*, and *C*) to the first capillary vessel segment (*order 0*) with harmonic wave frequency of 0.001 Hz and 1.5 Hz. Figure 11*B* shows the phase angle data for the respective coronary subtrees at 1.5 Hz.

In Fig. 12, the pressure and flow waves were calculated along the trunk of LAD artery from the inlet to the first capillary with simulated blood viscosity [based on the viscosity model of Pries et al. (26)]. Figure 12*A* shows the flow waves sequentially at 2.0-cm intervals along the trunk starting from the inlet of LAD artery. Figure 12*B* shows the pressure waves sequentially at different diameters along the trunk of LAD artery.

To examine the effect of viscosity, we maintained the same inlet and outlet boundary conditions and considered a uniform viscosity (μ = 1.1 cp) a compared with a blood viscosity distribution that depends on diameter (26). The wave frequency was maintained as 1.5 Hz in both cases. Figure 13 shows the relationship between mean segment flow |Q|_{mean} (average of the flow rate at every order) and mean segment pressure |P|_{mean} {average of 0.5 × [|P(0,ω)| + |P(*L*,ω)|] at every order} for the two viscosity models. As shown in Fig. 13, the flow rate becomes smaller at the same pressure when the blood viscosity is increased.

## DISCUSSION

#### Low-frequency flow (ω → 0) compared with previous steady-state flow model.

For any long cylindrical vessel, the input impedance for zero frequency can be found as follows: (11) where *L*(*n*) and *D*(*n*) are the length and diameter of the vessel segment, respectively; ω is the angular frequency; *Z*(0,0) and *Z*(*L*,0) are the impedance at the entrance and exit of the vessel, respectively; and μ is blood viscosity. It is known that the Womersley number, α = (*D*/2), describes the ratio of transient inertia force to viscous force. When α(1) << 1, the viscous force dominates and the flow can be assumed to be steady-state flow. In the present computation of the low-frequency flow (ω → 0), the Womersley number α(1) ≈ 0.07, 0.055, and 0.07 for the inlet vessel segment of the LAD, LCX, and RCA trees, respectively. The Womersley number at the inlet vessel segment is greatest for the coronary arterial tree. *Equation 11* and the Womersley's analysis suggest that the pulsatile flow should be steady-state flow when the frequency of the pulsatile flow approaches zero. Figure 4 shows that the pulsatile flow results agree reasonably well with the steady-state except that there is a small deviation below *order 4* because of slightly different outlet boundary conditions between the steady-state and pulsatile models. This comparison suggests that the current mathematical model predicts the correct results for the special case of steady-state flow.

#### Experimental validation of mathematical model.

To validate the present mathematical model, the numerical simulation must be compared with experimental data. Figure 5 shows the distribution of experimentally measured mean pressure and flow ± SD at the inlet of LCCA and RCA of six hearts. In Figure 5, it is observed that the mean flow rate at the inlet of LCCA is about twice as large as that of RCA, which is in agreement with previous studies (30). In Fig. 6, the predicted flow rate by the mathematical model agrees reasonably well with experimental results (within ± 1 SD). Therefore, the present mathematical model correctly simulates the pulsatile flow in diastole in the absence of vessel tone in the complex anatomical tree, which has more than 1 million vessels. It should be noted that the experimental preparation contains both the arterial and venous trees, whereas the current simulation only considers the arterial tree. It appears that the arterial tree dominates the wave form, however, given the present agreement. This may be due to the small diameter of the capillary network, which strongly retards the propagation of the flow wave. Furthermore, the resistance of the arterial tree is significantly larger than that of the venous counterpart.

#### Comparison with other studies.

The present mathematical model is compared with Olufsen's model for the root impedance (24, 25). In the present model, the constitutive equation is derived from the experimental data for the coronary arteries (see details in appendix), which may be written as follows: (12) where *E*(*n*)*h*(*n*)/*R*(*n*) is obtained directly from the linear curve fit of the experimental data. The dynamic elastic properties of the coronary artery are considered in the formulation (see details in appendix). The constitutive equation in the Olufsen's model (24, 25) is written as: (13) where *E*(*n*)*h*(*n*)/*R*(*n*) = *k*_{1}exp(*k*_{2}*r*_{0}) + *k*_{3} (see Refs. 24, 25). Figure 14 shows a comparison of the measured pressure-cross-sectional area (pressure-CSA) relationship between these two models. Figure 7 shows a comparison of the flow waves when these two models are implemented in the entire coronary arterial trees (22). Although there is a phase angle and amplitude difference between the results, it is found the mean flow rates are similar. The reason for the phase angle difference is due to the difference in *E*(*n*)*h*(*n*)/*R*(*n*).

Figure 7*B* also shows the results of the two models (present and Olufsen; Refs. 24, 25) compared with Zamir's formulation when extended to the entire tree (6, 33, 34). A significant difference (a factor of 5) appears between Zamir's analysis and those of Olufsen and the present study. In Zamir's mathematical formulation, the reflection coefficient *R*(*k,j*) at the downstream junction of bifurcation has the form (*Eq. 27* of Ref. 6): (14) where *Y*(*k,j*) is the characteristic admittance of the mother vessel segment (*k,j*). *Y*_{e}(*k* + 1,2*j* − 1) and *Y*_{e}(*k* + 1,2*j*) are effective admittances at the entrance of the two daughter vessels [(*k* + 1,2*j* − 1) and (*k* + 1,2*j*)], respectively. The difference stems from the application of the characteristic admittance *Y*(*k,j*) rather than the effective admittance at the exit of the mother vessel. Although this error may be negligible for small networks (such as in Ref. 6), it is significantly propagated in a network of millions of vessels. Hence, the error may be accumulated when all effective admittances (see Ref. 6) are calculated by marching backward along the tree structure and hence the observed discrepancy. It should be noted that the Zamir (6) and Olufsen (24, 25) formulations were not intended for the present vascular tree structure or the present capillary boundary conditions, and hence a direct comparison is not entirely possible.

#### Effect of various parameters on pulsatile blood flow.

In pulsatile flow, the wave frequency may have an effect. The pressure and flow waves at different heart rates produced by the piston pump are thus analyzed at the inlet of the LAD artery. In Fig. 8, the measured and computed flow wave show similar waveform at different heart rates. To ensure accurate comparison, the mean pressures (as shown in Table 1) are fixed at the same value. The computed mean flow rates are very similar to those of steady-state flow because the zero frequency in Fourier series has a dominant effect.

Because the study of harmonics is important, the dependence on frequency of flow wave is examined. In Table 2, the real part of the complex flow rate [Q_{inlet}(0,ω)] does not change with the wave frequency even if the flow wave is harmonic with the same inlet pressure [P_{inlet}(0,ω)]. The phase angle of the flow wave is increased with increase in wave frequency. The amplitude |*Z*_{inlet}(0,ω)| and the phase angle of the impedance *Z*_{inlet}(0,ω) depict similar tendency as those in Ref. 24, as shown in Fig. 9. Naturally, the amplitudes |*Z*_{inlet}(0,ω)| are different because our anatomical coronary arterial trees are different from those used in Ref. 24. On the basis of the above analysis, we conclude that the wave frequency has a small effect on the mean flow rate and has a significant effect on the phase angle and the flow waveform.

The significant asymmetry of coronary arteries is important to understand the spatial heterogeneity in coronary flow. At present, the flow wave along the trunk and major branches of the coronary arterial tree is analyzed with the harmonic wave frequencies of 0.001 Hz (∼steady state) and 1.5 Hz (90 beats/min). It is noted in Fig. 10, *B* and *A*, that the blood flow through the trunk shows a gradual drop along the path to the capillary vessels while the flow from the trunk to the primary branches has an abrupt drop because the trunk vessel segments have much larger diameter than the primary branches. The difference of the amplitudes |Q(*x*,ω)| is not apparent between the harmonic wave frequencies of 0.001 Hz and 1.5 Hz. However, the phase angle of Q(*x*,ω) at the harmonic wave frequency of 1.5 Hz reaches 36° at the inlet of the trunk and gradually decreases along the trunk. The phase angle of Q(*x*,ω) of the primary branches (i.e., *branches A*, *B*, and *C*) shows a nearly linear decrease toward the capillaries.

Most of the previous studies on pressure wave transmission have concentrated on changes along the aorta and into the lower limb arteries. This is the first investigation of the pressure and flow wave transmission throughout the entire coronary arterial tree. Because the patterns of waves in the LAD, LCX, and RCA trees are similar, the pressure and flow waves are only shown for the trunk of LAD artery from the inlet to the first capillary. As seen in Fig. 12*A*, changes in amplitude of the flow wave are apparent at 2.0-cm intervals along the trunk because the primary branches shunt the flow away from the trunk. A small phase angle shift of the flow waves is gradually developed along the trunk from the inlet to the first capillary. Figure 10*B* also predicts that the pressure only has an abrupt drop when the diameter is below 100 μm, which is in agreement with experimental measurements (2). Figure 12*B* depicts the relation between pressure waves and diameter along the trunk of LAD artery. It is noted that the attenuation of the amplitude of the pressure wave is rather small for larger vessels and becomes significant for vessels <100 μm. The phase angle shift of the pressure wave shows a similar tendency to that of the flow wave.

As shown in Fig. 13, the viscosity has an important effect on the blood flow in the entire coronary arterial trees, particularly in the lower-order vessels (smaller diameter). The viscosity retards the movement of blood flow. The Womersley number may be written as α = (*D*/2)() = transient inertia force/viscous force. An increase in the viscosity of pulsatile blood flow will increase the effect of viscous force so that the flow rate will decrease as predicted by the model.

#### Critique of model.

Although our model is based on detailed measured morphometric data of coronary vessels, there are still a number of assumptions made. For example, *1*) the ad hoc resistance boundary conditions are imposed on capillary vessels to decouple the arterial from the venous system; *2*) all vessels in the present model are assumed to be in vasodilated state in the absence of vessel tone where coronary flow reserve is substantially reduced; *3*) in the absence of experimental data, we assume that the measured compliance of larger vessels (19) applies to the entire vasculature; and *4*) the Womersley-type model is an approximate solution and does not incorporate nonlinear advective losses.

In future studies, the present arterial circuit will need to be extended to the entire coronary vasculature; i.e., the capillaries need to be connected to the entire venous system. Hence, the ad hoc resistance boundary condition on the capillary vessel will need to be replaced with a measured boundary condition at the outlet of the venous system. The present mathematical model must ultimately include the myocardial/vascular interaction during systole. Future studies are also needed to determine the mechanical properties of small coronary arterial vessels in situ. On the computational side, the one-dimensional/three-dimensional equations must be solved by using the impedance boundary condition based on the anatomical trees.

#### Significance of model.

The present study yields a mathematical model of pulsatile blood flow that incorporates some of the detailed anatomical features that influence coronary blood flow. The accuracy of the model is determined through experimental validation. The validated model of normal hearts will serve as a physiological reference state. Pathological states such as hypertension and stenosis can then be studied in relation to changes in model parameters. The proposed study makes use of physical principles, with the help of anatomy and experimental measurement, to explain and predict some aspects of pulsatile nature of the coronary arterial circulation in quantitative terms. The present model will serve to quantitatively test various hypotheses in the coronary circulation under pulsatile flow conditions.

## APPENDIX

In the present study, the blood flow is assumed to be impermeable, incompressible, Newtonian, and laminar. The blood flow in a vessel segment is taken as the flow in an axisymmetric cylinder. Therefore, the velocity of the fluid inside the vessel is denoted by *u*_{r}(*r,x,t*) and *u*_{x}(*r,x,t*), where *r* is the radial coordinate, *x* as the position along the length of vessel, *t* is time, *u*_{r} is the radial velocity, and *u*_{x} is the axial velocity. The cross-sectional area of the vessel *A*(*n*) is calculated from the diameter of the vessel *D*(*n*) (*n* is the vessel number) consistent with experimental measurements. The density ρ and viscosity μ are considered to be constant.

Wave propagation in a tube (e.g., the wave in a vessel *segment AB* as shown in Fig. 1) is governed by wave equations for the pressure p(*x,t*) and volume rate of flow q(*x,t*). To derive the wave equations, the continuity (mass conversation) and momentum equations for tube flow may be written as follows: (A1) (A2) It is noted that there are two equations and three independent variables [*A*(*x,t*), q(*x,t*) = 2π∫_{0}^{R(n)}*u*_{x}(*r,x,t*)*r*d*r*, and p(*x,t*)]. Therefore, an additional equation in the form of constitutive relation is needed to solve for the three unknown variables.

The constitutive equation is based on a pressure-CSA relation for every vessel. From Laplace's law, we have: (A3) where τ_{θ} is the mean circumferential stress that relates the transmural pressure p(*x,t*) − p_{0} (intravascular pressure minus external pressure, which is assumed zero in diastole), the vessel radius *R*(*n*) = *D*(*n*)/2, and the wall thickness *h*(*n*) for each order of vessel *n*. If we assume that the stress-strain relation is linear, we may write: (A4) where *E*(*n*) is Young's modulus for each order vessel and [*R*^{2} − *R*^{2}(*n*)]/*R*^{2}(*n*) is the circumferential strain. If we combine *Eqs. A3* and *A4*, we obtain the following constitutive equation: (A5) where the radius is expressed in terms of cross-sectional area. At this point, a pressure-CSA relation is needed. Such relations were previously measured by our group for the coronary arteries proximal to 1.0 mm in diameter using digital subtraction angiography (19). It was found that the differences in the pressure-CSA relation of vessels proximal to 1.5 mm in diameter (*orders 11, 10, 9*) were relatively small (see Fig. 3 and Table 1 in Ref. 19). Therefore, the mean of the pressure-CSA data for the 1.3 mm and 2.8 mm diameter vessels (see Fig. 3 in Ref. 19) were adopted here. As determined by linear least-squares fits of data, the experimental data may be expressed as follows: (A6) which can be approximated as: (A7) If we combine *Eqs. A5* and *A7*, we obtain *E*(*n*)*h*(*n*)/*R*(*n*) = 4.0 × 10^{5}. The thickness-to-radius ratio for the 1.3- and 2.8-mm vessels is very similar and was used to calculate the Young's modulus *E*. When *R*(*n*) and *h*(*n*) are obtained from Guo and Kassab's study (13), the static Young's modulus *E*(*n*) is determined as ∼8.0 × 10^{6} (dyn/cm^{2}). In the present simulations, it was assumed that this value is constant through the coronary arterial tree.

The dynamic Young's modulus must be considered because the current model calculates the impedance, pressure, and flow under various frequencies after a Fourier transform. Previous studies (4, 12) have investigated the dynamic elastic properties of canine coronary artery. In the current study, the ratio of *E*_{dyn}/*E*_{stat} (4) is used to adjust the dependence of Young's modulus on frequency.

By assuming the blood flow is harmonic and quasi-steady state, the variables can be written as: (A8) (A9) (A10) (A11) where ω is the angular frequency. Also, P(*x*,ω), *U*_{x}(*r,x*,ω), Q_{x}(*x*,ω), and *A*(*x*,ω) may be simplified as P, *U*_{x}, Q, and *A*. If we substitute *Eqs. A5–A11* into *Eqs. A1* and *A2*, we obtain: (A12) and (A13) The velocity profile was first obtained by Womersley (31) for pulsatile flow in a rigid tube as a solution to *Eq. A13* in the form: (A14) where α = *R*(*n*) is the Womersley number and *J*_{0} is the zero-order Bessel function. The volume flow rate is the same as that obtained by Duan and Zamir (5). Integrating the flow velocity over the cross-sectional area, we obtain the following wave equation: (A15) where *F*_{10}(α) = [2*J*_{1}(*i*^{3/2}α)]/[*i*^{3/2}α*J*_{0}(*i*^{3/2}α)] and *J*_{1} is the first-order Bessel function. Finally, by using *Eqs. A12* and *A15*, the wave equations for the pressure P and volume rate of flow Q may be written as follows: (A16) (A17) where *c*_{0} = is the wave velocity without the viscous effect, and *Y*_{0} = *A*(*n*)/ρ*c*_{0} is the characteristic admittance. Also we define *Z*_{0} = 1/*Y*_{0} as the characteristic impedance. If we combine Eqs. *A16* and *A17*, we obtain: (A18) When the viscous effect is incorporated into the wave velocity, we may define the wave velocity *c* = ·*c*_{0}. *Equation A18* can then be written as follows: (A19) Defining *Y*_{1} = *Y*_{0} and *Z*_{1} = *Z*_{0} and solving *Eqs. A16* and *A19* yield *Eqs. 1* and 2 and subsequently *Eqs. 3–7* as outlined in mathematical model.

## GRANTS

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

## Acknowledgments

We thank Dr. X. Lu and Dr. Y. Huang for assisting in the experimental preparation.

## 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 © 2006 by the American Physiological Society