Heart and Circulatory Physiology

A hybrid one-dimensional/Womersley model of pulsatile blood flow in the entire coronary arterial tree

Yunlong Huo, Ghassan S. Kassab


Using a frequency-domain Womersley-type model, we previously simulated pulsatile blood flow throughout the coronary arterial tree. Although this model represents a good approximation for the smaller vessels, it does not take into account the nonlinear convective energy losses in larger vessels. Here, using Womersley's theory, we present a hybrid model that considers the nonlinear effects for the larger epicardial arteries while simulating the distal vessels (down to the 1st capillary segments) with the use of Womersley's Theory. The main trunk and primary branches were discretized and modeled with one-dimensional Navier-Stokes equations, while the smaller-diameter vessels were treated as Womersley-type vessels. Energy losses associated with vessel bifurcations were incorporated in the present analysis. The formulation enables prediction of impedance and pressure and pulsatile flow distribution throughout the entire coronary arterial tree down to the first capillary segments in the arrested, vasodilated state. We found that the nonlinear convective term is negligible and the loss of energy at a bifurcation is small in the larger epicardial vessels of an arrested heart. Furthermore, we found that the flow waves along the trunk or at the primary branches tend to scale (normalized with respect to their mean values) to a single curve, except for a small phase angle difference. Finally, the model predictions for the inlet pressure and flow waves are in excellent agreement with previously published experimental results. This hybrid one-dimensional/Womersley model is an efficient approach that captures the essence of the hemodynamics of a complex large-scale vascular network. The present model has numerous applications to understanding the dynamics of coronary circulation.

  • coronary flow
  • pulse wave
  • admittance
  • impedance

the fluid dynamic approach uses the geometry and mechanical properties of the vasculature and the principles of conservation of mass and momentum to obtain the blood flow-pressure relation. The solution of such equations yields information on instantaneous velocity and pressure distributions. The fluid dynamic approach was first formulated in 1755 by Euler. Womersley (32, 33) and McDonald [as reported by Nichols and O'Rourke (22)] linearized Euler's equation and presented the theory and solution in great detail. Because of the complexity of the vascular architecture, however, this approach recognizes that it is necessary to idealize the geometry of the vasculature.

The frequency-domain Womersley-type approach was recently applied to the entire coronary arterial tree down to the first segment of (arterial) capillaries in a vasodilated, potassium-arrested heart (7). Inasmuch as this approach represents a linearization of the fluid dynamic (time-domain) approach, it ignores higher-order nonlinear terms. Although this approximation may be reasonable for the smaller vessels, the nonlinear effects must be considered for the larger vessels, where the Reynolds and Womersley numbers are higher (14). Hence, for accurate prediction of the distribution of transient pressure and flow at any position in the main trunk and primary branches, the vessels should be discretized and the time-domain fluid dynamics method must be used.

There are two major time-domain approaches to simulate the pulsatile blood flow and pressure waves in the vascular system: the “one-dimensional” (1-D) model of Hughes and Lubliner (6) and the “wave-intensity analysis” model of Parker et al. (25, 27, 30, 31). Since the high computational cost becomes prohibitive if the entire coronary arterial tree is considered, the models require appropriate outflow boundary conditions to simulate a truncated tree. Different outflow boundary conditions have been developed to satisfy the time-domain equations in lumped distal models, such as the pure resistive load boundary conditions (29), the windkessel model (10, 17), the nonreflecting outlet boundary conditions (1), or the impedance boundary conditions (23, 24). Despite the utility of lumped models, they provide limited insight into regional issues (5). Hence, the objective of the present study is to develop a hybrid 1-D/Womersley model that not only incorporates the nonlinear effects in the larger epicardial vessels but, also, integrates the smaller vessels down to the arterial capillaries.

Here, the 1-D model of Hughes and Lubliner (6) is used for the study of pulsatile blood flow in the main trunk and primary branches, along with the impedance boundary conditions. The impedance at each junction between the large and small vessels (time vs. frequency domain) is calculated through Womersley-type numerical computations (7) on the basis of the morphometric data of the full coronary arterial tree (15). The implicit finite-difference solver is developed to solve the 1-D model. The predictions of the mathematical model are compared with the previous Womersley-type numerical analysis and the experimental measurements on vasodilated, potassium-arrested porcine hearts. The agreement is found to be very good. The limitations and physiological implications of this hybrid model are discussed.


Governing equations.

The details of the mathematical derivations are outlined in appendix a. Briefly, the governing equations for flow and pressure may be expressed through conservation of mass and momentum as Math(1) Math(2) where A is cross-sectional area (CSA) of the vessel, q is volumetric flow rate, ρ is density, Estat is static Young's modulus, h is wall thickness, R0, h0, and A0 are original radius of the vessel, original wall thickness, and original CSA, respectively, and ν is kinematic viscosity. As outlined in appendix a, Eq. 2 can be rewritten as follows Math(3) The constitutive equation for a coronary vessel can be expressed as (7) Math(4) where p and p0 are the internal and external pressures, respectively.

Equations 1, 3, and 4 are used to calculate the pulsatile blood flow in each segment of the larger arteries (e.g., the main trunk and primary branches; Fig. 1, A and B). Viscosity was considered constant at 1.1 cP to mimic the cardioplegic solution used in the experiments. We shall establish relevant boundary conditions to extend the model to the entire coronary arterial tree.

Fig. 1.

Schematic representation of computational domains (main trunk and primary branches) in right coronary artery (RCA, A) and left anterior descending coronary (LAD) and left circumflex (LCx) arterial trees (B). C: schematic of branching angles.

Boundary conditions.

The inlet pressure boundary condition is obtained from experimental measurements (7). Olufsen's impedance boundary conditions (23, 24) are adopted as the outlet boundary conditions at the terminals of the 1-D model (Fig. 1, A and B). Womersley's theory from a recent publication (7) is applied to the morphometric trees to represent the distal vascular beds for each of the outlets of the numerical domain. Impedance/admittance [Z(x,ω)/Y(x,ω)] is calculated at each outlet of the numerical domain. By inverse Fourier transformation, z(x,t)/y(x,t) may be obtained from Z(x,ω)/Y(x,ω). With use of the convolution theorem, the new outflow boundary conditions may be obtained as follows Math or Math(5) which is used to prescribe the pressure-flow boundary conditions at each outlet. For the junction boundary condition, mass is conserved at each junction Vortices created at the bifurcations can result in loss of energy. When the effect of gravity is neglected, a loss coefficient K (18, 24, 27) can be incorporated into Bernoulli's equation as follows Math(6) Math(7) where ux is axial velocity. Equations 6 and 7 can be used to determine the pressure-flow relation at each junction in the large coronary arteries.

Branching angles.

To calculate the loss coefficient, the optimum branching angles (Fig. 1C) are calculated using formulations described elsewhere (21, 34). The optimal branching angles are given as a function of area ratios as Math and Math(8) where Math(9) with R0 BC < R0 BD. R0 BC and R0 BD are the radii for vessel segments BC and BD (Fig. 1C), respectively. Once the optimum branching angles are calculated, the loss coefficients can be estimated from Ref. 18.

Material parameters.

The viscosity (μ) and density (ρ) of the solution were selected as 1.1 cP and 1 g/cm3, respectively, to mimic our experimental cardioplegic solution containing albumin. The coronary wall thickness for every order was adopted from data of Guo and Kassab (3). The static Young's modulus was calculated as ∼7.0 × 106 dyn/cm2, as described in a previous publication (7).

Method of solution.

To solve the nonlinear hyperbolic 1-D equations, we adopted the time-centered implicit (trapezoidal) finite-difference method, which is stable and second order in time and space. The details of the mathematical derivations are outlined in appendix b. Since the method is stable, we checked only the effect of different meshes in the present study. The criterion for the convergence of the iteration between different meshes was set to 1 × 10−3 (norm-2 relative error). Finally, the mesh size (Δx) was selected as 0.05 cm, and the time step (Δt) was set to 2.0 × 10−3. The run time for the FORTRAN program is ∼1 h when the general LU decomposition method is used to solve the sparse matrix assembled in the largest computational domain on an AMD Opteron 240 computer. Even for a large sparse matrix with millions of equations, the matrix can be solved relatively quickly by using the LU decomposition with partial pivoting and triangular system solvers through forward and backward substitution (SuperLU_dist implemented in ANSI C and MPI for communications; Ref. 9).


The computational model described above was employed to predict the transient pressure and flow waves in the main trunk and primary branches. Extensive numerical simulations have been carried out for various parameters and computational domains. A selection of computed results is compared with the previous experimental and theoretical results (7). The computational domains in the 1-D model are shown in Fig. 1, A and B, for the right coronary artery (RCA), left anterior descending coronary artery (LAD), and left circumflex artery (LCx). The 1-D model is restricted to the trunk of the major arteries (proximal artery to ∼2-mm-diameter vessel).

Pressure boundary conditions at the inlet of the LAD and LCx.

It is desirable to verify the time-domain computer code for the values of pressure and flow with available experimental or numerical solutions. In the present study, the flow waves calculated using the hybrid 1-D/Womersley model with pressure and impedance boundary conditions at the inlet and outlet of the RCA, LAD, and LCx are compared with the experimental measurements and the full Womersley model (7). Figure 2, A and B, shows the flow wave at the inlet of the RCA, LAD, and LCx, as shown in Fig. 1, A and B, respectively. Comparison with experimental results shows slightly better agreement with hybrid model results than with Womersley model results. Both models are within ±1 SD of the experimental data (Fig. 2).

Fig. 2.

Pulsatile flow in hybrid 1-dimensional (1-D)/Womersley model and full Womersley-type model and experimental results at the inlet of RCA (A) and LAD and LCx (B). Both models are implemented in the entire coronary arterial tree.

The 1-D model was also simulated without the second-order viscous term. The resulting flow waves (data not shown) were indistinguishable from those in Fig. 2, which include the viscous term. The computed time-averaged ratios of pressure to velocity head and convective to pressure term at the inlet of vessel segments are listed in Table 1 for various branch angles. The pressure term is much larger than the convective term. Furthermore, the ratio of convective to pressure term is much smaller for the primary branches than for the trunk vessels.

View this table:
Table 1.

Various parameters for proximal LAD vessel segments in Fig. 1B

The hybrid model can accurately predict the transient blood flow in the main trunk and primary branches. The flow waves were calculated along the main trunk of the LAD. Figure 3 shows the flow waves sequentially at different spatial positions along the main trunk starting from the inlet of the LAD. The decrease in the amplitude is apparent, but the flow waveform remains relatively unchanged, except for a very small phase angle shift in the potassium-arrested heart. Figure 4 shows the flow waves at the inlet of various primary branches. When the data are normalized relative to the mean of each waveform, the patterns appear remarkably similar at every primary branch. Figure 5, A and B, shows the flow waves normalized by the mean value at the inlet of the main trunk and primary branches, respectively. Figure 5C shows the normalized flow waves at the inlet of the main trunk and one primary branch. The patterns of the normalized flow waves are very similar, except for a phase angle shift. The change of the amplitude of the pressure waves in the 1-D computational domain is small, because attenuation of the amplitude of the pressure waves is significant for <100-μm vessels. Therefore, we do not show the pressure waves along the main trunk.

Fig. 3.

Flow waves sequentially along the main trunk (TA–TB) starting from the inlet of the LAD.

Fig. 4.

Flow waves of primary branches (PA, PB, and PC) of the LAD trunk (TA–TB).

Fig. 5.

A: flow waves normalized by time-averaged value in main trunk of LAD. B: flow waves normalized by mean value at inlet of primary branches of LAD. C: normalized flow waves at inlet of main trunk and 1 primary branch of LAD.

Figure 6 shows the relation between mean stem flow (average over time) and stem diameter in the main trunk and primary branches of the RCA, LAD, and LCx.

Fig. 6.

Relation of time-averaged stem flow to stem diameter in main trunk and primary branches of the RCA, LAD, and LCx. Solid line corresponds to least-square fit of data according to a power law relation, with correlation coefficients of 0.98, 0.97, and 0.97 for RCA, LAD, and LCx, respectively, and exponents of 2.40, 2.30, and 2.13, respectively.

Effect of energy loss at a bifurcation.

Since vortices or separation of flow often occur at bifurcations, the energy loss should be considered. The loss of energy at a daughter vessel segment of a bifurcation is affected by many parameters, such as the branching angles (θ), the Reynolds number (Re), the area ratios (A/Amother), and the flow rate ratios (Q/Qmother). Table 1 lists these various parameters for the vessel segments in the first seven bifurcations. Various loss coefficients were estimated in the range 0.25–0.45, which yields almost identical flow waves compared with zero loss coefficient. Figure 7 shows a comparison of flow waves when loss coefficients are 0, 0.35, and 1.5. The computed flow wave shows a slight apparent dissipation as the loss coefficient increases. Obviously, the dissipation increases when a larger portion (beyond the 7 bifurcations) of the tree is considered.

Fig. 7.

Transient flow when the first 7 bifurcations are considered. K, loss coefficient.


A novel hybrid 1-D/Womersley model.

In a recent study (7), pulsatile blood pressure and flow distribution in the entire coronary arterial tree down to the first capillary segments was analyzed using a full Womersley-type model. Although the Womersley-type model reveals the spatial and temporal flow and pressure distribution in the entire coronary arterial tree, it neglects many important aspects, such as nonlinear convective losses and the diffusion term in Navier-Stokes equations. In particular, it is difficult to obtain a dynamic solution during cardiac contraction from a Womersley-type model to evaluate the myocardial interaction.

In the present study, we focus on solving the 1-D model of Hughes and Lubliner (6) in large epicardial vessels (the main trunk and primary branches; Fig. 1, A and B) in the potassium-arrested heart. The impedance boundary conditions are computed through a Womersley-type model at every flow outlet. Both models are based on the measured morphometric data of Kassab et al. (15, 20). This is the first investigation to implement the time-domain model coupled with the frequency-domain model in the entire coronary arterial tree on the basis of experimentally measured morphometric data of the entire coronary arterial tree.

The present agreement between the hybrid and Womersley models implies that the nonlinear convective term plays a minor role, even in the larger epicardial vessels in an arrested heart. Table 1 shows that the pressure term is ≥100 times larger than the convective term in the Navier-Stokes equation and that the effect of the pressure head outweighs the effect of the velocity head in the Bernoulli equation. We also found that the viscous term is small, because solution of the 1-D model with and without the second-order viscous term yielded very similar results. A stable solution is reached after several iterations after assignment of initial values of every parameter (e.g., the initial flow and initial CSA for every vessel).

The pressure and flow waves have been widely studied from the aorta to the limb arteries (22). For the first time, the pressure and flow waves were calculated along the main trunk of the LAD in the present study. As described recently (7), the amplitude of the flow waves in Fig. 3 decreases gradually, and a small phase angle shift of the flow waves appears along the main trunk of the LAD in the potassium-arrested heart, because the primary branches shunt the flow away from the main trunk. Since primary branches have different CSA, the flow waves at the inlet of primary branches show different amplitudes. When the flow waves are normalized by the time-averaged flow rate, however, the flow waves tend to scale to a single curve, except for a small phrase angle difference. This novel observation reflects an interesting scaling of flow to structure.

The observation that the normalized flow waves maintain similar patterns relates to the scaling laws of coronary circulation (13). Specifically, we found that Murray's law holds for the present analysis, but with an exponent of 2.4, 2.3, and 2.1 for the RCA, LAD, and LCx, respectively. This value for the pulsatile analysis is slightly larger than that of previous steady-state analysis (2.1) primarily because of the higher viscosity used in the previous simulation. Future studies will attempt to generalize the steady-state Zhou-Kassab-Molloi model (35) to pulsatile conditions.

Large asymmetric bifurcations show a predilection for atherosclerosis. Such bifurcations have been extensively studied biologically (2, 4, 11, 12, 16, 26). Here, we include the loss coefficients to consider the loss of energy at a bifurcation in the 1-D model. The loss coefficients were estimated from Ref. 18 through the various computed parameters on the basis of the morphometric data. The loss coefficients are larger in the daughter vessels with smaller diameters, because those vessels have larger bifurcation angles than mother vessels. The loss coefficients were also applied to multiple bifurcations in the coronary arteries. The flow waves show small dissipation over time when the loss coefficients are increased in the 1-D model (Fig. 7). Hence, we found that the loss of energy at the bifurcation is small in the diastolic coronary arteries.

Critique of the model.

Although our hybrid 1-D model is derived from the Navier-Stokes equations with fewer assumptions than the Womersley-type model and is based on the detailed morphometric data of coronary vessels, a number of issues deserve mention. The time-domain impedance has large variation, which may cause instability when the 1-D computational domain is expanded. Conversely, small ripples or oscillation may appear if the computational domain is too small (restricted to a few bifurcations). This may be due to the impedance boundary condition. The time-domain impedance has a large variation when the inverse Fourier transformation is used to transform the frequency-domain Womersley-type impedance from the anatomic data to the time domain. A small error in the impedance boundary condition at any outlet of the thousands of vessel segments connecting to the large epicardial artery may cause instability. To remedy this difficulty, the imposed impedance boundary conditions should be validated with experimental measurements. To further improve the accuracy, the second-order trapezoidal finite-difference method should be implemented with a variable time increment option and higher-order Runge-Kutta formulation.

Significance of the model.

The present time-domain model is capable of accurately predicting the dynamic pulsatile flow on the basis of detailed coronary morphometric data. The 1-D model can be made patient specific through interface with patient-specific anatomy of large epicardial vessels obtained from CT, MRI, or other imaging techniques. The smaller vessels that are beyond clinical imaging resolution can be treated as lumped. This approach will significantly contribute to the design of diagnostic, therapeutic, and treatment devices, such as catheters, grafts, and stents. When mass transport is considered, the model can further predict the distribution of drugs, cells, and genes to design treatment for the ischemic myocardium.

Future studies.

In the future, the present mathematical model will be implemented in a 3-D coronary arterial tree (11) and coupled with a 3-D myocardial mechanical model to determine the myocardial-vascular interaction during systole. Specifically, the hybrid model will predict the dynamic pressure and flow distribution in the larger coronary arteries during systole when the time-dependent impedance or admittance is calculated through the Womersley-type model coupled with regional stresses and strains in the contracting myocardium. Furthermore, 2-D and 3-D numerical models will be developed to depict the detailed flow streamlines. The 1-D model will be adopted as a dynamic boundary condition, instead of the lumped parameters, for future 2-D/3-D simulation of flow streamlines. This will allow the prediction of flow at the critical sites, where vortices or separations may occur. Furthermore, the shear stress and more accurate loss coefficients induced by vortices or separation of flow will be quantified. Finally, the pressure-induced vibration of the large epicardial vessel wall will be considered in 2-D and 3-D models to solve the solid-fluid structure through the arbitrary Lagrangian-Eulerian formulation (front-tracking technique) (1, 8).


Derivation of 1-D fluid dynamic equations for large epicardial arteries.

The 1-D equations are formulated for a Newtonian, incompressible fluid in an elastic tree. The 1-D analysis consists only of larger coronary arteries, every vessel of which is assumed to be cylindrical, with an impermeable wall. Therefore, the fluid flow in every vessel segment is axisymmetric and laminar with a no-slip boundary condition (i.e., the velocity of fluid at the wall equals the velocity of the wall). Various geometrical parameters (transient/initial physiological measurements) in a typical vessel can be represented as length (L/L0), radius (R/R0), CSA (A/A0), surface (S/S0), volume (V/V0), and wall thickness (h/h0). Furthermore, the dynamic parameters are represented as the velocity of fluid flow [ur(r,x,t), ux(r,x,t)]. The volumetric flow rate and the pressure in the vessel are q(x,t) and p(x,t) (intravascular pressure − external pressure, which is assumed zero in diastole), respectively. A cylindrical coordinate system was used with radial (r) and axial (x) directions; the hemodynamic quantities ur(r,x,t), ux(r,x,t), q(x,t), and p(x,t) may be represented as ur, ux, q, and p.

To derive the 1-D equations, the continuity (mass conservation) and momentum equations for tube flow (6) may be written as follows Math(A1) Math(A2) where ν = μ/ρ is kinematic viscosity. Since we assumed the no-slip boundary condition and deformation of the elastic wall of the vessel, we obtain the following boundary equations Math(A3) The definition of volumetric flow rate and CSA is Math(A4) Math(A5) Integrating Eqs. A1 and A2 over the CSA of the vessel with the boundary conditions (Eq. A3), with the assumption that ρ is constant over the area, we obtain Math(A6) Math(A7) In the larger epicardial arteries, the velocity profile is assumed to be fully developed Poiseuille flow. We assume the parabolic velocity profile, thus obtaining the following equation Math(A8) where C1, C2, and C3 are integration constants. When the boundary conditions (Eq. A3) are considered, we solve Eq. A8 to obtain Math(A9) Using Eqs. A9 and A4, we find Math(A10) The second term of Eq. A7 may be written as Math(A11) and the viscous drag force (right side of Eq. A7) can be expressed as Math(A12) Inserting Eqs. A10 and A11 into Eq. A7, we obtain Math(A13) From our previous study (7), the constitutive equation has the form Math(A14) where Estat is obtained from the experimental results. Inserting Eq. A14 into Eq. A13, we obtain Math(A15) Subsequently, Eqs. 1 and 2 are outlined in mathematical model.


Method of solution.

Equations 1 and 3 can be summarized as follows Math(B1) Math where Math For the large vessels, Eq. B1 is solved by the time-centered implicit (trapezoidal) finite-difference method. We define q Math= q(mΔx,nΔt) and AMath= A(mΔx,nΔt), where 0 < nN for the current time step and 0 < mM expresses the position along a vessel divided into M subintervals. Therefore, the discretization in Eq. B1 may be written as follows Math(B2) Math(B3) where Math Math

Inlet boundary condition.

The inflow q into the coronary arteries (RCA, LAD, and LCX) is obtained from experimental measurements. The inflow q Math and qMath can be obtained directly from the measurements. When the pressure pMath and pMath are given, the real area AMath and AMath at the inlet of the coronary arteries can be calculated on the basis of Eq. 4. If inflow or inlet pressure is given, the other may be calculated. For example, if the inflow qMath and qMath are given, we obtain the following equations at the inlet Math(B4) Math(B5) Math(B6) When the inlet pressure p Math and pMath is given (e.g., in the cases of the LAD and LCx), we obtain the following equations at the inlet Math(B7) Math(B8) Math(B9)

Outlet boundary condition.

The convolution integral (Eq. 5) can be written as follows Math or Math(B10) which can be discretized as Math or Math(B11) At the outlet, we obtain the following equations Math(B12) Math(B13)

Junction boundary condition.

The simplest and most prevalent junction is a bifurcation (15), which represents an outflow boundary for the mother vessel and an inflow boundary condition for the daughter vessels (d1 and d2). For a bifurcation, at time n + 1, Eqs. 6 and 7 may be discretized as follows Math(B14) and Math(B15) The procedure can be summarized as follows Math(B16) Math(B17) Math(B18) Math(B19) Math(B20) Math(B21) Equations B16 and B17 are the continuity and momentum equations, respectively, at the outlet of the mother vessel segment; Eqs. B18 and B19 are the continuity and momentum equations, respectively, at the inlet of one daughter vessel segment, while Eqs. B20 and B21 are the continuity and momentum equations, respectively, at the inlet of the other daughter vessel segment. Equations B2B21, which constitute the formulation to solve the 1-D model, are in turn used at every bifurcation in the main trunk and primary branches.


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


We thank Henry Yu Chen for the three-dimensional Pro/Engineer CAD models of the LAD, LCx, and RCA.


  • 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.


View Abstract