Math-Tech and Department of Mathematics, Roskilde University,
Roskilde 4000, Denmark
A central
problem in modeling blood flow and pressure in the larger systemic
arteries is determining a physiologically based boundary condition so
that the arterial tree can be truncated after a few generations. We
have used a structured tree attached to the terminal branches of the
truncated arterial tree in which the root impedance is estimated using
a semianalytical approach based on a linearization of the viscous
axisymmetric Navier-Stokes equations. This provides a dynamic boundary
condition that maintains the phase lag between blood flow and pressure
as well as the high-frequency oscillations present in the impedance
spectra. Furthermore, it accommodates the wave propagation effects for
the entire systemic arterial tree. The result is a model that is
physiologically adequate as well as computationally feasible. For
validation, we have compared the structured tree model with a pure
resistance and a windkessel model as well as with measured data.
 |
INTRODUCTION |
THE ABILITY TO PREDICT blood flow and pressure at any
position along the larger systemic arteries can lead to a
better understanding of the arterial function. Wave shapes of arterial
pressure are strongly influenced by the wave reflections resulting from
tapering of the vessels, bifurcations, changes in the arterial
distensibility, and a high resistance in the arterioles, i.e., a high
downstream peripheral resistance (27, 31, 36). The appearance of the dicrotic wave can be important in clinical situations. For example, the
dicrotic wave is diminished in some patients suffering from diabetes or
vascular diseases such as atherosclerosis (13, 22, 23, 27, 31).
Furthermore, it has been observed that patients with stiffer arteries
have a less pronounced dicrotic wave but an increased systolic pressure
(5, 18, 27, 31). Therefore, studies of the dicrotic wave and comparison
of pressure profiles at different positions could possibly be used for
diagnostic purposes, e.g., to locate stenosis (39).
The purpose of this study was to formulate a nonlinear physiological
model predicting blood flow and pressure at any position along the
larger systemic arteries. These quantities are regarded as functions of
time and one spatial dimension; hence, by definition, the model is one
dimensional. We focus on how to restrict the computational domain while
maintaining a physiological approach, i.e., on deriving a
physiologically correct downstream boundary condition allowing the
model to be implemented and executed without excessive computational
costs. In previous arterial models (see, e.g., Refs. 1, 6, 12, 40,
45-47, 53) the peripheral beds have received little attention.
Models of the peripheral beds have mostly used rather simple boundary
conditions, e.g., pure resistor conditions (1, 47) or a windkessel
condition (38, 40, 45, 46).
The problem with these models is that they are lumped models, and thus
they cannot include wave propagation effects in the part of the
arterial system that they model. However, the overall wave shape can be
approximated by having good values for the total resistance and
compliance, but these quantities are not easy to measure and the pulse
profiles are sensitive to the value of these parameters. Furthermore,
it is hard to avoid artificial reflections in a system using lumped
models as outflow boundary conditions.
Our approach to overcoming these problems is to terminate the larger
arteries by attaching a structured tree representing the remainder of
the arterial system. This is illustrated in Fig. 1. Although the model for the larger
arteries is fully nonlinear, we construct a semianalytical solution
based on a linear hydrodynamic model for the structured tree.

View larger version (21K):
[in this window]
[in a new window]
|
Fig. 1.
Systemic arterial tree. Tree consisting of larger arteries, in which
nonlinear equations are solved, originates at the heart (A) and
terminates at the points marked with B. Structured trees representing
smaller arteries originate at these terminals and provide main tree
with outflow boundary conditions. r, Radius;
ri, i = 0-3 is the minimum
radius used for terminating the structured trees.
|
|
Physiologically it makes sense to split the model in two parts. The
role of the larger arteries is to distribute blood to all parts of the
body, whereas the role of the smaller arteries is to allow perfusion of
specific tissues. In fact, several papers (see, e.g., Refs. 9, 44)
showed that the smaller arteries are distributed in an optimal and
structured way such that they cover the tissue evenly using a
minimization principle. They also showed that the larger arteries do
not follow such rules. Furthermore, blood flow in the larger systemic
arteries is dominated by inertia, whereas blood flow in the smaller
arteries is dominated by viscosity (11).
This paper is divided into four sections, the model for the larger
systemic arteries and its boundary conditions, the model for the
smaller arteries (the structured tree model), the results of these
models, and discussion. It should be emphasized that relations derived
for the larger arteries do not necessarily apply to the smaller
arteries because of the inherent differences mentioned above.
 |
LARGER SYSTEMIC ARTERIES |
A one-dimensional model for blood flow and pressure in the larger
systemic arteries can be derived using Navier-Stokes equations. Numerous papers discuss such models (see, e.g., Refs. 1, 7, 40, 46,
47). They differ in the way they treat the shear stresses, the relation
between pressure and cross-sectional area, and the boundary conditions.
The derivation presented here is based primarily on Barnard et al. (7)
and Peskin (36), but where appropriate we discuss some of the other
approaches. The larger systemic arteries can be viewed as compliant and
tapering vessels connected in a bifurcating tree. We have chosen to
define the larger arteries as those shown in Fig. 1. For computational reasons we lump some of the branches together (the coronary arteries, the intercostals, and the arteries branching from the celiac axis), and
we assume more symmetry, by letting arteries that have both a left and
right branch be identical, than is actually the case in the human body.
We assume that the fluid is incompressible and Newtonian, that the flow
is axisymmetric and laminar, and that the vessels are circular and
tapering, i.e., that they can be represented by a surface of
revolution (see Fig. 2) with length
L, radius R, cross-sectional area A, surface S,
and volume V. Although this holds for the larger arteries (14), it
does not apply to the smaller arteries and arterioles because they do
not taper significantly. Furthermore, we assume that S
moves with velocity v = [vr(x, t),
vx(x, t)], the fluid
moves with velocity u = [ur(r, x, t),
ux(r, x, t)],
and the pressure in the system is
p(r, x, t).
R(x, t) is the actual radius of the vessel, and
r and x are the cylindrical coordinates in radial (r) and length (x) directions.

View larger version (13K):
[in this window]
[in a new window]
|
Fig. 2.
A typical vessel shown in a cylindrical coordinate system
(r, radial coordinate; x, longitudinal
coordinate; t, time). R(x,t) is the
radius, and A(x,t) is the
cross-sectional area. The end surfaces are given at x = 0 and x = L (vessel length).
|
|
From Barnard et al. (7) and Peskin (36) we get the following
axisymmetric Navier-Stokes equations for conservation of volume and
x momentum
|
(1)
|
|
(2)
|
where
is the density, µ is the dynamic viscosity, and
the last term represents the viscous wall shear stress. Equations 1 and 2 constitute the exact mathematical model for the
system. It is not one dimensional, because both
ux and p vary with r as well as
x. To obtain a one-dimensional model we assume that p and thus
p/
x are functions of x and t only, i.e.,
that the pressure is constant over the cross-sectional area.
Furthermore, we assume that the velocity profile is flat but with a
boundary layer of thickness
, hence
where
is the mean value of
ux(r,x,t).
According to Lighthill (24) the thickness of the boundary layer for the
larger arteries can be estimated from (
/
)1/2 = [
T/(2
)]2
0.01 cm, where
= µ/
= 0.046 cm2/s is the kinematic viscosity,
is the
angular frequency, and T = 1.25 s is the period of one
heartbeat. Measurements show that the velocity profile changes
throughout the arterial system, from being almost flat in the aorta to
a more parabolic shape in the peripheral arteries (25, 33, 34). This is
a result of the fact that the boundary layer remains constant (1 mm) as
the vessels are getting smaller.
Integrating over the cross-sectional area and disregarding terms of
first and higher orders of
gives
|
(3)
|
|
(4)
|
where q(x, t) = A(x, t)
(x, t)
is the flow and p(x, t) is the pressure averaged
over the cross-sectional area.
Equations 3 and 4 are the basic equations for the
one-dimensional model of arterial pulse wave propagation. However, we
only have two equations and three dependent variables, namely,
p(x, t), q(x, t), and
A(x, t). Therefore, we need a third
relation, the so-called state equation. This is based on the compliance of the vessels and gives an equation for the motion of the vessel walls. Several approaches can be taken on how to model these
"elastic" properties. The arterial wall is not purely elastic but
exhibits a viscoelastic behavior (11, 25, 42), i.e., there is a time delay in the response from a change in pressure to the corresponding change in cross-sectional area. However, to keep the model simple we
disregard viscoelasticity and use a relation derived from the linear
theory of elasticity. This is a reasonable assumption because the
viscoelastic effects are small within the physiological range of the
flow and pressure (50). We need to examine the equilibrium of the
internal and external forces acting on a unit element of the wall. We
assume that the arterial vessels are circular, that the walls are thin,
i.e., that wall thickness (h)/r
1, that the loading and deformation
are axisymmetric, and finally that the vessels are tethered in the
longitudinal direction. In this case the external forces are reduced to
stresses acting only in the circumferential direction (3) and from what
is often known as Laplace's law we get the circumferential tensile
stress
|
|
where pe = p
p0 is the excess
pressure, i.e., the pressure of the vessel minus the pressure of the
surroundings, h is the wall thickness,
(r
r0)/r0 is
the corresponding circumferential strain, E
is
Young's modulus in the circumferential direction, 
=
x = 0.5 are the Poisson ratios in the
circumferential and longitudinal directions, and r0
is the radius at zero transmural pressure, i.e., at p = p0.
Because we assumed that the vessel is tethered in the longitudinal
direction this is the only contribution we get when balancing the
internal and external forces, and we can without loss of generality
drop the
subscript on E. Solving for pe we get
|
(5)
|
where A0(x) =
r0(x)2 is the
cross-sectional area at zero transmural pressure. It should be noted
that Eq. 5 has the property that the area
A(x, t) becomes infinite at a finite
transmural pressure ("blow out"). Real arteries resist this
tendency by having a nonlinear Young's modulus E that
increases with increasing strain. For a given position x in
each arterial segment it is possible to compute
Eh/r0 from the corresponding radius
r0(x), estimates for the volume
compliance Cvol, and the approximation
Plotting Eh/r0 as a function of
r0 one can fit (empirically) the exponential
function
|
(6)
|
to these estimates. k1,
k2, and k3 are constants.
With data for Cvol from Westerhof et al. (53), Stergiopulos
et al. (46), and Segers et al. (45) we obtain k1 = 2.00 × 107
g · s
2 · cm
1,
k2 =
22.53 cm
1, and
k3 = 8.65 × 105
g · s
2 · cm
1. The
data and the fitted curve are shown in Fig.
3.

View larger version (17K):
[in this window]
[in a new window]
|
Fig. 3.
Young's modulus E times wall thickness (h) relative to
radius at zero transmural pressure (r0) as a
function of r0.
|
|
Boundary Conditions
The system of Eqs. 3-5 is hyperbolic, and the square of
the wave propagation velocity
(c20) is
Because c0 is always positive, and the
wave propagation velocity generally is much larger than the velocity of
blood
(|u|
c0),
the characteristics will cross and have opposite directions. We thus
need one boundary condition at each end of the vessels. The boundary
condition at the inflow to the arterial tree is applied at A in Fig. 1,
and the outflow boundary conditions are applied at each of the
terminals, marked by B in the figure. Finally, we need three conditions
at each of the bifurcations to close the system of equations.
Inflow boundary condition.
At the inflow boundary we need to specify the flow, the pressure, or a
relation between them. Because the shape of the pulse wave in the upper
ascending aorta is generated by the inflow from the aortic valve, we
have chosen to represent the inflow by a periodic extension of a
measured flow wave (see Fig. 4). This curve
is measured in the upper ascending aorta using magnetic resonance (20).
The measured curve is modified in a number of ways. First, it is made
periodic, such that q(0, 0) = q(0, T ) where T = 1.25 s is the length of the period. Second, it is smoothed to avoid too
many oscillations in our simulated results. Finally, we have scaled the
curve so that the cardiac output is changed from 4.14 to 4.03 l/min, a
reduction of 3%. This is a consequence of the fact that we have not
included all branches in our arterial model, and hence there is a
certain amount of blood not included in our system.
Outflow boundary condition.
The outflow boundary condition can be determined in several ways. The
simplest reasonable approach is to let the outflow be proportional to
the pressure, i.e., to let the boundary condition be determined by a
pure resistive load. This is commonly used in previous work (1, 29, 43,
47, 48). However, it is not obvious how to choose the correct value for
the peripheral resistance at the points where the larger arteries are
terminated. Furthermore, if we assume a constant relation between flow
and pressure at the downstream boundary they are forced to be in phase, which is generally not physiologically valid in these relatively large
arteries. This is also pointed out by Anliker et al. (2), who note that
the pure resistance boundary condition only applies if the arteries are
sufficiently small. This can be seen (see Fig. 7A) from the
hysteresis curves appearing when we plot p(x, t)
versus q(x, t) parametrized by t and for a
fixed x. The forced in-phase condition propagates back through
the vessel, changing the overall slope as well as narrowing the width
of the loop, which means that the phase is disturbed throughout the
vessel. Because we are looking for some reflections in the system (just enough to produce the dicrotic wave), we would expect a small change in
the loop rating curves but not as drastic as the one appearing with
this boundary condition.
As mentioned in the introductory paragraphs, the other approach often
used is to apply a windkessel model at the outflow boundary (see, e.g.,
Refs. 40, 46, 52), among others. The most commonly used windkessel
model is the three-element model, which represents the resistance and
elasticity of the vessels by an electrical analog model consisting of a
resistance in series with a parallel combination of a resistance and a
capacitor (see Fig. 5). The frequency-dependent impedance (Z ) of the windkessel model
is given by
|
(7)
|
where RT = R1 + R2 is the total resistance and CT is
the total compliance of the vascular bed. These three parameters must be specified for each boundary condition. Also from this model the
hysteresis curves, appearing when plotting p(x, t)
versus q(x, t) parametrized by t and for a
fixed x (see Fig. 7B) show that flow and pressure are
almost in phase. However, the windkessel model does not change the
slope and width of the curves significantly. Furthermore, none of these
models is able to include wave-propagation effects. Therefore, we have
investigated how the physical domain extends beyond the boundary of the
larger arteries using a one-dimensional fluid dynamic approach.

View larger version (7K):
[in this window]
[in a new window]
|
Fig. 5.
Three-element windkessel model. For each terminal vessel, where the
model is applied, 3 parameters must be estimated, that is, resistances
R1 and R2 as well as total
compliance (capacitance) CT.
|
|
The arterial system consists essentially of a large asymmetric tree
with a varying number of generations, ranging to >20 generations before the precapillary arterioles are reached. It would be too comprehensive to compute the full nonlinear model for such a tree. Therefore, a more appropriate strategy is to describe the flow and
pressure in these smaller arteries using a simpler model that can be
solved analytically, e.g., a linear model. From these subtrees for the
smaller arteries we can obtain a suitable boundary condition for the
system of nonlinear equations as a time-dependent relation between flow
and pressure. We treat these subtrees of smaller arteries as a
structured network of straight vessels in which the corresponding
linear equations are solved; this is treated in detail in SMALLER
SYSTEMIC ARTERIES. From these solutions it is possible, using
Fourier analysis, to determine a more physiological relation between
flow and pressure. In this case we actually see a phase lag between
flow and pressure (see Fig. 7C).
Because the inlet boundary condition is periodic, we assume that the
flow and pressure can be expressed using complex periodic Fourier
series. Then any feature of the system response can be determined
separately for each term. Let
where
k = 2
k/T is
the angular frequency and
For any Fourier mode we define the frequency-dependent
impedance Z(x,
) as
|
(8)
|
where we have used the terminology of electrical networks,
with P playing the role of voltage and Q the role of current. If we can
find an expression for Z(x,
) and hence by
inverse Fourier transform find z(x, t), we
arrive at an analytic relation between p(x, t) and
q(x, t) by the convolution theorem
|
(9)
|
This is our new outflow boundary condition for the larger
arteries, which should be evaluated at each of the terminals marked by
B in Fig. 1, i.e., at x = Li where
Li is the length of the ith
terminal segment. In SMALLER SYSTEMIC ARTERIES we discuss
in detail how the impedance can be obtained.
The convolution integral in Eq. 9 spans one entire period;
hence, it requires knowledge of the solution at all times. But, because
we are solving Eq. 3 using an explicit method (Richtmeyer's 2-step version of Lax-Wendroff's method), the solution will not be
computed simultaneously at all times in the period. However, because
the pulse wave propagation is a periodic phenomenon, this difficulty
can be overcome by using an iterative approach to determine the
convolution integral; at any time t we use values from the previous period for t'
[t:T ] and
the ones already calculated for t'
[0:t]. We
then iterate over a number of periods until a stable solution is
reached. Numerical experiments suggest that this happens after no more
than three to four periods.
Bifurcation conditions.
The bifurcations can be modeled using a number of different approaches
(1, 24, 47). We assume that the bifurcation takes place at a point;
hence, we need three conditions to close the system of equations.
Because no blood is leaving the system at that point we get
|
(10)
|
The subscript pa refers to the parent vessel, and the
subscripts d1 and d2 refer to the daughter
vessels. Assuming that the pressure is continuous over the bifurcation
|
(11)
|
we get the other two conditions. These conditions pose some
questions because of the Bernoulli law that may or may not apply depending on the details of the flow pattern at the junction. At a
boundary where the total cross-sectional area decreases proceeding downstream, one would, according to the Bernoulli law, expect a drop in
pressure associated with the increase in velocity. In the arterial
system, however, the total cross-sectional area typically increases at
junctions (again, proceeding downstream toward the periphery), and
hence an increase in pressure would be expected. On the other hand,
because the change in area at the junction is discontinuous, flow
separation and vortex formation are expected just downstream from the
bifurcation and the Bernoulli law does not apply. Under these
circumstances, which invoke dissipation of kinetic energy, it is more
appropriate to use pressure continuity.
 |
SMALLER SYSTEMIC ARTERIES |
Modeling biological networks as structured trees has been done
previously but mainly for the pulmonary airways or for smaller systems
of arteries, such as the coronary arteries (8). However, because the
purpose of all of the smaller arteries is to cover some tissue evenly
with blood we have found it natural to assume that the smaller arteries
as a whole are structured in some way.
We assume that it is possible to construct an asymmetric binary
structured tree where at each bifurcation the radius of the two
daughter vessels is scaled by factors
and
(0 <
,
< 1), respectively, and where the branching
terminates when the radius of any given vessel is less than some given
minimum radius. We then need to derive a system of equations that gives
the root impedance of the structured tree. Such a tree is shown in Fig. 6 and is indicated at the outflow from each
of the larger arteries in Fig. 1.

View larger version (17K):
[in this window]
[in a new window]
|
Fig. 6.
A structured tree. At each bifurcation radius of the daughter vessels
are scaled by a factors and , respectively, i.e., for any radius
rpa of a parent vessel the radii of its 2 daughters
are rpa and rpa. Because
each branch is terminated when the radius is less than some given
minimum radius, the tree does not have a fixed number of generations.
|
|
In the following sections we first set up the system of equations
needed to determine the root impedance and then discuss the actual
geometrical properties of the tree.
Root Impedance
We have chosen to use the linear hydrodynamic model originally
suggested by Womersley (54), Atabek and Lew (4), and Pedley (35). We
assume that blood is incompressible and Newtonian and that the
viscosity terms, which are dominant for the smaller arteries (11), are
included in the derivation of the flow equations. The model predicts
flow Q(x,
) and pressure P(x,
) in the
frequency domain by balancing the forces of the elastic wall with those acting inside the fluid. As for the larger arteries, the fluid equations are based on the axisymmetric Navier-Stokes equations for
flow in a circular cylinder and the wall equations are determined from
balancing forces on a thin shell representing an infinitesimal element
of the surface of the cylinder. The boundary condition linking these
equations together is the no-slip condition that adheres the fluid
particles to the inner surface of the tube and hence to the motion of
the elastic wall. The result after averaging over the cross-sectional
area and linearizing the equations is the following x momentum
equation
|
(12)
|
where
and w = r(
/
)1/2 is the
Womersley parameter [for details see Atabek (3) and Pedley (35)]. It
should be noted that K is not continuous. The reason for this
is that it is derived from asymptotic expansions for w
0 and w
; however, because the jump is small it is
possible to construct a smooth function for K for all
w. We define the compliance C by the approximation
The latter approximation applies because Eh
pr0 and the corresponding continuity equation
therefore becomes
|
(13)
|
Equations 12 and 13 can be combined into the
following wave equation
|
(14)
|
A similar equation can be derived for P. Solving Eq. 14 for Q and using this solution together with Eq. 12 or
13 gives P. From these we get the following equation for the
impedance
|
(15)
|
where a and b are arbitrary constants.
Furthermore, we have introduced the wave propagation velocity c =
(which is complex) and the short-hand notation g =
. If we
assume that Z(L,
) is known, it is possible to
find an expression for b/a using Eq. 15
Evaluating Eq. 15 at x = 0 and inserting the
expression for b/a gives the following input impedance
|
(16)
|
To extend this analysis to the more general case of a tree
we need to determine how to "cross" a bifurcation.
This does not require much more than we have already established. As
for the larger arteries, we assumed in Eqs. 10 and 11
that the pressures are continuous and that the flows, and hence the admittances, add
|
(17)
|
If we know g and L for each vessel together
with the impedance at each of the distal terminals (leaves of the tree)
we can use Eqs. 16 and 17 recursively to find the input impedance.
The peripheral resistance for the larger arteries arises mainly from
viscous flow in the arterioles (27). In our model of the smaller
arteries (the structured tree) we include sufficiently many branches on
the arteriolar level to generate a realistic peripheral resistance.
Hence we can assume a zero impedance at the leaves of the structured
tree. However, arterioles are muscular and are able to dynamically
regulate the flow to the organs in question depending on their need.
This can be modeled by changing the zero impedance to a nonzero
impedance at the leaves of the structured tree. Furthermore, it is
possible to simulate vasodilation or vasoconstriction by changing the
radius and the elasticity of the vessels. Vasodilation, for example,
could be simulated by lowering k1 in the relation
for Eh/r0 (6), implying a lower value for
the very small radius, as well as increasing the radius at all levels
of the structured tree but without changing the asymmetry ratio or the
length of the vessels. This can be done for some or just one of the
vessels because the parameters for the relation are local to each branch.
For any vessel, the average impedance, or in electrical terminology the
DC term (for
= 0), can be found as
|
(18)
|
where
is a constant defining the length-to-radius ratio;
this is discussed further in Geometry of Structured Tree. Equation 18 suggests that in general for any network the root impedance will
be proportional to r
30. For the case we are looking at, where the tree is terminated whenever the
radius of the leaves of the tree is less than some given minimum radius, the constant of proportionality cannot be derived analytically, but in the special case of a symmetric tree with N generations we get
Geometry of Structured Tree
To construct the structured tree in which Eqs. 16 and 17 are solved for all the individual vessels, we need
information about the order of the tree and the impedance at the
terminals and some structured way to present the geometry of the
vessels, i.e., the length, the radius, the wall thickness, and Young's modulus.
From knowledge about the total cross-sectional area of the arterial bed
we can estimate the order of the structured tree. The arteriolar
diameters of the vascular bed vary from 10 to 50 µm according to the
organ in question. At a diameter of 50 µm (for large arterioles) the
total cross-sectional area is ~55 cm2 (27). In case of a
symmetric tree this gives ~2.8 × 106 terminal branches.
Because we start approximately at the second generation after the aorta
and create a number of structured trees this must be modified
accordingly. In case of the tree presented in Fig. 1 there are 21 subtrees, and we therefore need ~17 generations, because
However, the arterial tree is asymmetric, and therefore this
estimate will only apply approximately. Therefore, the best condition
to apply is to continue the binary bifurcation of the structured tree
until the diameter of any vessel in the structured tree becomes less
than the arteriolar diameter of the vascular bed in the given organ.
At each bifurcation we need a law determining how the geometry (radius
or cross-sectional area) changes over the junction. Such a relation is
suggested by Uylings (51) among others. It is derived from the
principle of minimum work in the arterial system and is given by
|
(19)
|
where rpa is the radius of the parent
vessel, and the subscripts d1 and d2 refer to
the two daughter vessels. This relation is valid for a range of flows;
the radius exponent
= 3.0 is optimal for laminar flow and
= 2.33 for turbulent flow. The exponent has been discussed rather
extensively in the literature (37, 49), and on the basis of these
reports we have chosen
= 2.76.
We could choose to construct a completely symmetric tree. Although this
is not what appears in the physical world, it may still reflect the
essential behavior of the tree. In a symmetric tree all vessels will be
terminated at the same point, and as a result the impedance propagated
from the leaves of the tree will be in phase accentuating the
reflections from the boundary. This will not happen in an asymmetric
tree. In this case the impedance is scattered, and as a result
reflections will be attenuated. Hence, it makes sense that the arteries
in the human body branch asymmetrically. Furthermore, we can easily
include an asymmetry relation in the structured tree, i.e., the ratio
between the cross-sectional area of two daughter branches. The
following relations for the area and asymmetry ratios are suggested by
Zamir (55)
|
(20)
|
where
and
are constants that characterize the
structured tree. The parameters
,
, and
are not independent
but are related as follows
|
(21)
|
As with the radius exponent
, the values for
have
also been discussed extensively (see especially Refs. 16, 32, 51). On
basis of these reports we have chosen to keep
constant throughout the tree at a value of 1.16, giving an asymmetry ratio of
= 0.41. Using these relations we are able to determine the structure of the
tree, i.e., how the radius of the two daughter vessels should scale
relative to their parent. Assuming that they scale with factors
and
, respectively, these can be found from Eqs. 19 and 20 as
It is also necessary to describe the length of each vessel
as a function of some property of its parent, because this parameter appears in Eq. 16. There are a number of papers discussing how to determine the lengths of the vessels, but only a few of these suggest any connection with the other geometric parameters. Iberall (17) estimates that the length-to-radius ratio
lrr = l/r
50 ± 10. This is also supported by Suwa et al. (49), Bergel (10), and
Bassingthwaighte and Van Beek (9), among others.
Finally, to compute the compliance C and wave speed c, we need
information about how wall thickness h and Young's modulus E change or depend on other parameters throughout the tree.
This is exactly what we have estimated for the larger arteries, and because the smaller arteries are essentially composed of the same types
of tissue, we have extended the dependencies to apply here as well,
i.e., Eh/r0 is estimated as a function of
r0 as shown in Fig. 3.
Together these considerations constitute all the information needed to
construct an asymmetric binary structured tree. The tree will be
structured with respect to the radius in such a way that all parameters
are related to the radius of a given vessel and at each bifurcation the
radius of the daughter vessels scale with
and
,
respectively, where 0 <
,
< 1.
 |
RESULTS |
The results in this section are generated from the model shown in Fig.
1, in which we solve Eqs. 3-5 in the larger arteries, with
the inflow condition from Fig. 4 and the outflow condition (8) obtained
by solving Eqs. 16 and 17 in the structured tree. These
are then compared with both the windkessel and the pure resistance
models as well as with experimental data.
First, we present some features from the outflow boundary condition,
and afterwards we give some results showing the overall behavior of
this model.
Outflow Boundary Condition
To judge the performance of the structured tree model we will show
1) a comparison of the three outflow boundary conditions (the
pure resistance condition, the windkessel model, and the structured
tree model) applied to a single isolated vessel and 2) a
comparison of measured data for the impedance in humans and results
from simulations using both the windkessel and the structured tree
boundary conditions. The advantage of both the pure resistance and
windkessel models is that they are easy to understand and computationally inexpensive, whereas the disadvantage is that they are
not able to capture the wave propagation phenomena in the part of the
arterial system that they model. Furthermore, neither the pure
resistance nor the windkessel model can account for the phase lag
between flow and pressure. The windkessel model requires estimates of
the total arterial resistance, RT = R1 + R2, and compliance
CT for each terminal segment, and the pure resistance model
needs the total arterial resistance. Still, when coupled to the
nonlinear equations for the larger arteries, both models are able to
capture the overall behavior of the system.
To show the differences between the three models we have (for
simplicity) used a single tapering vessel of length 100 cm, with top
radius 0.4 cm and bottom radius of 0.25. We then have applied the three
outflow boundary conditions to the vessel. To make them match as well
as possible we have estimated the parameters for the windkessel and the
pure resistance models from the root impedance determined by the
structured tree model. RT (the DC term from the
structured tree model) is the same for all three models, and
CT for the windkessel model is fitted empirically to match
that given by the structured tree model. It should be emphasized that
this study is theoretical, and hence the parameters should not be
compared with physiological values. However, the same differences can
be seen when applying the three models to the whole tree.
The result of this comparison is shown in Fig.
7, in which we have plotted pressure versus
flow at five equidistant locations along the vessel. Figure 7A
is for the pure resistance model, Fig. 7B for the windkessel
model, and Fig. 7C for the structured tree model. When these
figures are compared, the most striking difference is that the pure
resistance model affects the overall shape of the curve. The forced
in-phase condition at the outflow boundary results in a narrowing of
the width of the loop back through the vessel. Furthermore, it is worth
noticing that the pure resistance model can be seen as a special case
of the windkessel model incorporating only the DC resistance. For the
windkessel model flow and pressure are also nearly in phase, but the
narrowing is not reflected back through the vessel. Finally, it is
observed that the structured tree model does indeed retain some phase
lag between flow and pressure. The overall shape of the pressure
profiles are similar, but there are some significant differences. These have to do with the fact that the structured tree model includes wave
propagation effects for the entire tree, which the windkessel model
cannot do. Hence, the windkessel and pure resistance models are likely
to introduce artificial reflections. When we compare the result for the
pressure waves (see Fig. 8), this is
exactly what we see as a difference between the models. The reflections from the pressure when the windkessel model and especially the pure
resistance model are used are more pronounced than for the structured
tree model.

View larger version (24K):
[in this window]
[in a new window]
|
Fig. 7.
Pressure (p; mmHg) vs. flow (q; cm3/s) over 1 cardiac cycle
for a single vessel at 5 equidistant locations, i.e., for a vessel of
length L, with x = 0, L/4, L/2,
3L/4, and L, where curve with highest flow is for
x = 0. Both q and p are plotted as functions of time
during 1 cardiac cycle. A: pure resistance model. B:
windkessel model. C: structured tree model.
|
|

View larger version (39K):
[in this window]
[in a new window]
|
Fig. 8.
Pressure in a single vessel with 3 outflow boundary conditions as a
function of x and t during 1 period. A: pure
resistance condition. B: windkessel condition. C:
structured tree condition.
|
|
However, these differences are more easily seen if we instead make a
so-called Bode plot of the impedance at the bottom of the large vessel
[ZL(
)] versus the angular frequency (
),
i.e., making a plot of
log(|ZL(
)|) versus
log (
) and a plot of phase (Z ) versus log(
). In this
paper we have made such plots for the windkessel and the structured
tree model and compared these with human measured data [adapted from
Nichols and O'Rourke (27)]. However, all these are from the larger
arteries. Therefore, to be able to compare we have applied the two
models directly as outflow boundary conditions for these larger
arteries even though the outflow boundary conditions usually are
applied further downstream. We show the comparisons of the windkessel,
the structured tree, and the measured data for the brachiocephalic
artery. However, similar results can be obtained for the other parts of
the arterial tree (both larger and smaller arteries). It should be
noted that the structured tree model is not designed to be valid for
the large arteries, and as a result one should not assume perfect matches without some adjustments of the parameters. For the
brachiocephalic artery we had to modify the length-to-radius ratio from
50 to 130. This much larger length-to-radius ratio corresponds with the
arteries of the arm. From anatomic data (see, e.g., Ref. 27), we see
that starting from the subclavian artery no large side branches occur
before the bifurcation between the ulnar and interosseous arteries,
which, in fact, gives a length-to-radius ratio of ~130. The
brachiocephalic artery has a major bifurcation after only ~3.5 cm,
resulting in a very short length-to-radius ratio. However, this is then
followed by a rather long length-to-radius ratio that is not taken into
account here. Because we know that the length-to-radius ratio gets
smaller for the smaller arteries, further studies might be able to
reveal some functional dependence, e.g., on the radius, of the
length-to-radius ratio. For the model of the brachiocephalic artery we
have chosen a root radius of 0.5 cm and a minimum radius of 0.025 cm. A
root radius of 0.5 cm is rather small, but Nichols and O'Rourke (27)
do not specify exactly where the data have been measured. Thus the
comparisons should be interpreted with all of these reservations in
mind. For the vessel wall parameters we have used the relation
Eh/r0 = k1
exp (k2r0) + k3 as shown in Fig. 3. The results for these comparisons are shown in Fig. 9.

View larger version (20K):
[in this window]
[in a new window]
|
Fig. 9.
A: brachiocephalic impedance
|Z(0, )| vs. angular frequency = 2 f for [0:125] Hz. B and C:
Bode plots, i.e., a log-log plot of modulus of the impedance
|Z(0, )| versus angular frequency = 2 f (B) and a single-log plot of phase of
impedance vs. frequency (C). Dotted lines indicate results
from windkessel model (wk), solid lines indicate results from
structured tree (st), and dashed lines indicate measured results (md 0 and md 1). Measured data are from Mills et al. (26).
|
|
From these results the difference between the models becomes evident,
namely, that the windkessel model cannot include the high-frequency
oscillations present in real human data. These observations can also be
seen when investigating the pressure profiles (see Fig. 8) in which the
reflections and the maximum pressure are more damped for the structured
tree model than for the windkessel or constant resistance models.
It is common knowledge that organs have different peripheral
resistances reflecting the physiological characteristics of the organs
they supply, and for any given organ or muscle the minimum radius
should be chosen to match the total peripheral resistance of the
corresponding organ. For example, the renals have a small peripheral
resistance reflected in a large minimum terminal radius, and the
femoral arteries have a large peripheral resistance reflected in a
small minimal terminal radius. For the simulations in this paper we
have used four different minimum radius values. These are
r0 = 0.06 cm, r1 = 0.04 cm,
r2 = 0.02 cm, and r3 = 0.01 cm. They are distributed as shown in Fig. 1. The compliance of the vessels
is determined from the pressure cross-sectional area relation using the
radius-dependent relation for Eh/r0.
Similar relations are needed for the windkessel model; it requires
estimates of arterial resistance and compliance.
This leads to the next question, namely, what happens if the root
radius is changed. This is interesting because the various structured
trees applied at the boundaries of the larger arteries have different
root radii. Because the most significant change appears in the average
impedance Z(0, 0), or in the electrical terminology the DC
term, we can assume that it should depend on the area, and
hence the radius, of the vessels in the structured tree. As expected
from our derivation in the impedance section the average impedance is
approximately inversely proportional to the root radius to the third
power. The agreement is not exact, however, because the number of
generations in the network is not fixed. Instead, the tree is continued
until a terminal radius is reached. Thus increasing
r0 results in a tree with more branches. The reason
for the decrease in impedance with increased root radius is that the
total cross-sectional area at the bottom of the structured tree is
increased, and hence it can accommodate a larger outflow. The total
cross-sectional area is increased simply because starting at a larger
root radius and terminating the tree when the radius of the terminals
are less than some fixed minimum value gives a larger tree with more terminals.
Coupled Model
Using the structured tree as a boundary condition is a feasible way of
determining the blood flow and pressure in the larger arteries. To
verify this we have chosen to depict pressures in the aorta and the
subclavian and brachial arteries, because these graphs show all the
significant phenomena and because they can be compared with measured
data (see Figs.
10-12).
Corresponding results for the other branches shown in Fig. 1 behave
similarly. The geometrical data for this model are taken from Segers et
al. (45) except that we have restricted the concept of larger arteries
to consist only of those that are at most one generation away from the
aorta, iliac, and femoral arteries (see Fig. 1).

View larger version (24K):
[in this window]
[in a new window]
|
Fig. 11.
Pressure p(xf, t) in aorta (for
xf = 0, 10, ... , 50 cm; A) and
subclavian and brachial arteries (for xf = 0,
10, ... , 40 cm; B).
|
|

View larger version (23K):
[in this window]
[in a new window]
|
Fig. 12.
Measured data for pressure p(x, t) in aorta
(A) and subclavian and brachial arteries (B) are
from Ref. 41.
|
|
First, we observe that the simulated pressure profiles exhibit all the
qualitative characteristics that must be present in an arterial model:
1) the maximum pressure increases away from the heart toward
the periphery; 2) the steepness of the incoming pressure
profile increases toward the periphery; 3) the reflected dicrotic wave separates from the incoming pressure wave and is more
prominent toward the periphery; and 4) the incoming wave of the
ascending aorta is round whereas the incoming wave, at more peripheral
locations, is narrower because it is accentuated by the reflected wave.
Aside from the features mentioned above, which are obviously also
present in the measured data, we observe a good quantitative
correspondence between our simulated results and the measured data. The
systolic and diastolic pressures are similar, and in the ascending
aorta there is a shoulder before the maximum of the incoming wave that
disappears at locations further distally. The latter is caused by the
large wave propagation velocity that results in a reflected wave
returning before the peak of the incoming wave. However, it should be
noted that the results of the simulation depend strongly on the inflow
profile from the aortic valve, cardiac output, the geometry of the
arteries, the minimum radius, and the relation for Young's modulus,
the wall thickness, and the vessel radius (6). Because these parameters are not quoted by Remington and Wood (41), who refer to the measured
curve as stemming from a normal subject, discrepancies are likely to
occur, especially because the pulse curve generally varies
considerably, depending on the parameters mentioned above, even among
normal subjects. Even though we do not know these specific details we
do see a great potential with a model such as the one described in this
paper, and we are currently working on setting up a number of
experiments ourselves to validate the model further.
 |
DISCUSSION |
This work has shown the feasibility of limiting the computational
domain. On the basis of physiological information alone we have
succeeded in deriving a boundary condition for the nonlinear model
predicting blood flow and pressure in the larger systemic arterial tree
that is able to retain the high-frequency oscillations present in the
impedance spectra. This applies at all levels, even at the peripheral
level where the boundary condition is applied. We also showed that
retaining the high-frequency oscillations is not possible using the
simpler pure resistance or windkessel models. Applying a structured
tree gives