## Abstract

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.

- arterial modeling
- mathematical modeling
- arterial outflow boundary condition

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.

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** = [*v*
_{r}(*x*, *t*),*v*
_{x}(*x*, *t*)], the fluid moves with velocity **u** = [*u*
_{r}(*r*, *x*, *t*),*u*
_{x}(*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.

From Barnard et al. (7) and Peskin (36) we get the following axisymmetric Navier-Stokes equations for conservation of volume and*x* momentum
Equation 1
Equation 2where ρ is the density, μ is the dynamic viscosity, and the last term represents the viscous wall shear stress. *Equations1
* and *
2
* constitute the exact mathematical model for the system. It is not one dimensional, because both*u*
_{x} 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*u*
_{x}(*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 cm^{2}/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
Equation 3
Equation 4where 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 p_{e} = p − p_{0} is the excess pressure, i.e., the pressure of the vessel minus the pressure of the surroundings, *h* is the wall thickness, (*r* − *r*
_{0})/*r*
_{0} 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 *r*
_{0}is the radius at zero transmural pressure, i.e., at p = p_{0}. 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 p_{e} we get
Equation 5where *A*
_{0}(*x*) = π*r*
_{0}(*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*/*r*
_{0} from the corresponding radius*r*
_{0}(*x*), estimates for the volume compliance C_{vol}, and the approximation
Plotting *Eh*/*r*
_{0} as a function of*r*
_{0} one can fit (empirically) the exponential function
Equation 6to these estimates. *k*
_{1},*k*
_{2}, and *k*
_{3} are constants.

With data for C_{vol} from Westerhof et al. (53), Stergiopulos et al. (46), and Segers et al. (45) we obtain *k*
_{1} = 2.00 × 10^{7}g ⋅ s^{−2} ⋅ cm^{−1},*k*
_{2} = −22.53 cm^{−1}, and*k*
_{3} = 8.65 × 10^{5}g ⋅ s^{−2} ⋅ cm^{−1}. The data and the fitted curve are shown in Fig.3.

### Boundary Conditions

The system of *Eqs. 3-5
* is hyperbolic, and the square of the wave propagation velocity (
) is
Because *c*
_{0} is always positive, and the wave propagation velocity generally is much larger than the velocity of blood (‖*u*‖ ≪ *c*
_{0}), 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. 7
*A*) 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
Equation 7where *R*
_{T} = *R*
_{1} +*R*
_{2} is the total resistance and C_{T} 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. 7
*B*) 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.

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. 7
*C*).

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
Equation 8where 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
Equation 9This 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* = *L*
_{i} where*L*
_{i} is the length of the *i*th 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
Equation 10The subscript pa refers to the parent vessel, and the subscripts d_{1} and d_{2} refer to the daughter vessels. Assuming that the pressure is continuous over the bifurcation
Equation 11we 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.

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
Equation 12where
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* ≫ p*r*
_{0} and the corresponding continuity equation therefore becomes
Equation 13
*Equations 12
* and *
13
* can be combined into the following wave equation
Equation 14A 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
Equation 15where *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
Equation 16To 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
Equation 17If 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 *k*
_{1} in the relation for *Eh*/*r*
_{0} (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
Equation 18where λ is a constant defining the length-to-radius ratio; this is discussed further in *Geometry of Structured Tree. Equation18
* suggests that in general for any network the root impedance will be proportional to
. 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 cm^{2} (27). In case of a symmetric tree this gives ∼2.8 × 10^{6} 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
Equation 19where *r*
_{pa} is the radius of the parent vessel, and the subscripts d_{1} and d_{2} 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)
Equation 20where η and γ are constants that characterize the structured tree. The parameters ξ, η, and γ are not independent but are related as follows
Equation 21As 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*l*
_{rr} = *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*/*r*
_{0} is estimated as a function of*r*
_{0} 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, *R*
_{T} =*R*
_{1} + *R*
_{2}, and compliance C_{T} 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. *R*
_{T} (the DC term from the structured tree model) is the same for all three models, and C_{T} 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 7
*A*is for the pure resistance model, Fig. 7
*B* for the windkessel model, and Fig. 7
*C* 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.

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 [*Z*
_{L}(ω)] versus the angular frequency (ω), i.e., making a plot of log(‖*Z*
_{L}(ω)‖) 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*/*r*
_{0} = *k*
_{1}exp (*k*
_{2}
*r*
_{0}) +*k*
_{3} as shown in Fig. 3. The results for these comparisons are shown in Fig. 9.

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*r*
_{0} = 0.06 cm, *r*
_{1} = 0.04 cm,*r*
_{2} = 0.02 cm, and *r*
_{3} = 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*/*r*
_{0}. 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*r*
_{0} 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).

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 a physiological boundary condition in which we only need to estimate the minimum order of the structured tree in accordance with the total peripheral resistance of the organ/muscle to which it leads. This is a result of the viscous treatment of blood flow in the smaller arteries in the structured tree. We thus get a peripheral resistance at the terminal of the larger arteries, i.e., the root impedance from the structured tree, entirely from the structured tree itself, and no terminal resistance needs to be applied at its leaves. This is consistent with the observations (15) that the smaller arteries generate the peripheral resistance and not the capillaries. However, to model pathological situations such as vasodilation or vasoconstriction it is possible to apply a nonzero terminal resistance, change the radius throughout the structured tree, or modify the relation including Young’s modulus if required. Furthermore, it should be noted that it is possible to estimate and compare the parameters of a pure resistance and a windkessel element boundary condition from the structured tree (30). The structured tree provides an alternative to the windkessel boundary condition, and its advantage is that it is based on the fluid mechanics of the arteries, hence it naturally includes the wave propagation effects, which cannot be modeled by a lumped model. The only parameter that must be estimated separately for each structured tree is the radius at the root of the tree.

When plotting the pressure of the larger arteries seen in Figs. 10 and11, we not only retain the phase lag, as seen in Fig. 7
*C*, and the high-frequency oscillations evident from the impedance plots, but we also get realistic results. However, it is necessary to reduce cardiac output with ∼3% to achieve these results. This is a natural consequence of the fact that we do not include all branches of the arteries, and hence we should lower the cardiac output accordingly. As we show in results our model reflects the most important phenomena of arterial function. The maximum pressure and the steepness of the incoming pressure waves are increased distally, and the dicrotic wave separates from the incoming pressure wave and is more prominent at more peripheral locations.

## Acknowledgments

The author thanks Prof. C. Peskin, Courant Institute of Mathematical Sciences, New York University, New York, NY, Prof. T. J. Pedley, Department of Applied Mathematics and Theoretical Physics, Cambridge University, Cambridge, UK, and colleagues at Math-Tech and the Department of Mathematics, Roskilde University, Denmark, for support and very fruitful discussions throughout this work.

## Footnotes

Address for reprint requests: M. S. Olufsen, Math-Tech, Admiralgade 22, DK-1066 Copenhagen, Denmark.

This work is a result of the author’s PhD study at Roskilde University, which was supported by the Danish Academy of Technical Sciences and by High Performance Parallel Computing, Software Engineering Applications Eureka Project 1063.

- Copyright © 1999 the American Physiological Society