Modeling blood flow in the arterial system

The investigation of arterial blood flow has always been an important topic in physiology. Large proportion of the cardiovascular diseases (arteriosclerosis, vasoconstriction) affects the arterial part of the circulatory system. Physiologists place huge efforts in the improvement of the diagnosis and the treatment of these diseases. During the last decades fluid mechanics has become a powerful tool in the analysis of arterial blood flow. In the current paper a numerical approach for the calculation of one dimensional unsteady blood flow in the arterial network is presented. The continuity and momentum equations are solved using the method of characteristics (MOC). Furthermore a viscoelastic material model is applied to describe the behavior of the arterial walls. The network model of the arterial system is set up using 45 viscoelastic branches. The heart is modeled using an unsteady volume flow rate boundary condition. Arterioles and capillaries are treated as linear resistance boundary conditions. The parame
ters of the viscoelastic material model were fine-tuned by comparing the results of the calculation with previously measured blood pressure curves. The resulting model is capable for simulating arteriosclerosis in an arbitrary part of the arterial network. The viscoelastic material model and the calculation method are presented in detail. Results of the calculation are presented and discussed.


Introduction
Simulating blood flow in the arterial network is a highly complex task.It consists of numerous vessel branches with different length, diameter and stiffness (compliance).There are also several bifurcations in the system.Furthermore the arterial vessel walls show a so called "bioviscoelastic" behavior (Monos [11]).This property has two important effects.First, the arterial wall deformation is a function of the transmural pressure and the time.Hence the "history" of the wall has an effect on the current state.Second, a hysteresis loop can be observed by loading and unloading the vessel wall.
In the last decades several one dimensional network models were introduced to describe blood flow mechanisms in the arteries.However most of these models were solved using lumped parameter methods like the impedance method (Avolio [3], Stergiopulos et al [15], Westerhof and Noordergraaf [19]).These approaches are limited in terms of modeling wave propagation.Loops in the network cannot be treated this way either.Distributed models were also used for modeling the arterial system.MOC is a popular mathematical method in this field (Streeter and Wylie [17], Halász [8]).Advantage of this method is that the non-linear terms of the momentum equation are not neglected like in the impedance method.Furthermore, the MOC ensures the information to travel with the pulse wave velocity.Anliker et al. [2] and Rockwell [14] used this method for modeling the arterial network.In both works a tapered tube was defined for this purpose.The effect of the side branches were modeled with a continuous leakage applied on the tube.This basic model was further developed by Stettler et al. [16].However their model consisted of only eight main arterial branches which is too few for accurate blood flow modeling.Wang and Parker [18] constructed an arterial network model containing 55 branches.However the pulse wave velocity was assumed to be constant.Also the viscoelasticity of the arterial vessel wall was neglected.
Other distributed models than the MOC were also used in research work.Azer and Peskin [4] used a two-step Lax-Wendroff method for their arterial network model.Their contribution is interesting as far as they take the Womersley theory into account.On the other hand the effect of the viscoelastic arterial vessel walls is neglected.Reymond et al [13] developed a method on the basis of the arterial network from Stergiopulos et al [15].A finite difference method was used for calculating the blood flow velocities and the pressures.
In the current paper the model network of the arterial system is set up using the schematic diagram of Avolio [3].The model was modified and simplified, the resulting network model has 45 arterial branches.A viscoelastic material model is introducedthe so called Stuart model is used for this purpose.To the best of our knowledge this material model was not yet used for arterial blood flow modeling.The describing equations of the Stuart model are solved together with the continuity equation and the one dimensional momentum equation using the MOC.The resulting numerical model accurately describes blood flow in the arterial system.The viscoelastic material model was validated in Licskó and Bárdossy [10].A sample venous network model was created and calculated in Bárdossy and Halász [5].The arterial model network presented in the current paper is further development of the previous publications.The human heart is modeled as a volume-flow rate boundary condition.Arterioles and capillaries are taken into account as linear resistances.The resulting blood pressure curves of the calculation show a fair agreement with physiological data.Therefore the model is capable for simulating arteriosclerosis in an arbitrary part of the arterial system, which is a new achievement in the research field.In the first calculation series a gradual closure of the right external iliac was investigated.The results are presented and analyzed from the physiological point of view.

Methods
In order to create the one dimensional model of the arterial network an appropriate numerical approach is needed.The Department of Hydrodynamic Systems at the University of Technology and Economics developed a software (the so called Transient Simulator) for calculating one dimensional unsteady fluid flow in pipe networks.Both lumped and distributed parameter models can be solved using this algorithm.The lumped parameter solver is used for modeling parts like pumps, valves, pressure-or mass flow sources.Using the distributed parameter solver the long elastic pipes are modeled which connect the lumped model parts.A one dimensional mesh is created along the elastic pipes, in every node pressure and velocity are calculated and stored.The continuity equation and the one dimensional momentum equation are defined for these inner nodes.The loss term of the momentum equation is valid for laminar flows.This system of partial differential equations is solved simultaneously using the MOC.The method is described in detail by Streeter and Wylie [17] and Alastruey et al. [1].The characteristic grid is defined on the basis of the spatial resolution (see Fig. 1).Pressure and velocity data are known in the inner nodes L and R at timestep t i-1 .Negative and positive characteristic curves are launched from these nodes.At the intersection of these nodes (P) pressure and velocity can be evaluated.
Boundary conditions are satisfied with only one characteristic curve (negative/positive) at the (upstream/downstream) end of the tube.Pressure or velocity data can be applied as boundary conditions.One main advantage of MOC is the accurate representation of pressure waves: during the calculation the information travels exactly with the wave-propagation velocity between the nodes of each pipe.
The Transient Simulator itself is mostly used for the calculation of unsteady flow in drinking water pipe networks.For instance pressure waves after an abrupt pump failure can be well estimated (Hős [9]).The software was also used for modeling waste water flow in open channels.In the Transient Simulator the Hooke law is applied for modeling the pipe walls.Blood vessel walls cannot be described accurately using the Hooke law since it neglects the time dependence of the deformation and there are no hysteresis loops by repeatedly loading and unloading the vessel wall.Therefore the next step was to extend the Transient Simulator with a viscoelastic material model.The Stuart model was selected for this purpose since it is used for modeling viscoelasticity in polymertechnics (Bodor and Vas [7]).The schematic diagram of the Stuart model is shown in Fig. 2. The Kelvin-Voigt element (ε 2 ) represents viscoelastic deformation while the elastic spring element (ε 1 ) models the elastic deformation component.
The describing equations of the Stuart model have to be combined with the governing equations of the fluid flow.The following equations have to be solved simultaneously: 1. Momentum equation: 2. Modified continuity equation: 3. Equations of the Stuart model: where D 0 and D are the initial and instantaneous diameters, p transmural pressure, v axial velocity, h elevation, δ 0 wall thickness, ν kinematic viscosity, g gravity and ρ density of the blood, E 1 , E 2 and η 2 material properties of the vessel.Wavepropagation speed (a), diameter and strain are calculated using the Stuart model: ) where ε 2e and ε e are the viscoelastic and total deformation from the previous timestep, p e is the transmural pressure from the previous timestep and t is the timestep.All of the above equations are solved within the Transient Simulator.All the applied approximations, assumptions and the detailed solution of these equations are discussed in Bárdossy and Halász [5].The model was validated with several measurements carried out on a single silicone tube (Licskó and Bárdossy [10]).After tuning the parameters of the Stuart model the agreement between measured and the simulated pressure curves was appropriate.
Using this model one dimensional unsteady flow in a network consisting of viscoelastic branches can be calculated (Bárdossy and Halász [5], [6]).The arterial network model was set up using these viscoelastic branches.The schematic diagram was constructed using the contribution of Avolio [3] and Wang and Parker [18].Some necessary modifications were applied and the network itself was simplified to 45 arterial branches.Discussions with physicians pointed out that this resolution is satisfactory for investigating arteriosclerosis in an arbitrary point of the arterial network model.The arterial branches are connected with their ending nodes to each other.In the bifurcations nodal continuity equations are defined and the static pressure is held constant.Inflow/outflow velocity of the connected branches can be calculated this way.The schematic diagram of the one dimensional arterial model is shown in Fig. 3 The heart was modeled using an unsteady volume flow rate source.The volume flow rate as a function of time was taken from Nichols and O'Rourke [12] and it is presented in Fig. 4.
The arterioles and the capillaries were modeled using resistances with linear behavior.
with R resistance and Q volume flow rate.Blood was modeled as a Newtonian fluid with an average density of 1050 kg/m 3 .
Setting proper initial conditions for each arterial branch and the interconnecting nodes is highly difficult.Therefore zero velocity, zero deformation and constant pressure was set.After starting the calculation the quasi-steady state is achieved within a few heart cycles.One main problem was setting up the parameters of the Stuart model for every arterial branch.In case of a single viscoelastic tube the tuning of the parameters E 1 , E 2 and η 2 can be carried out using a simple genetic algorithm.This method is able to fine tune the parameters of the Stuart model to fit measurement data.In case of the arterial network model three parameters have to be tuned for each branch.Since the network model consists of 45 branches, 135 parameters have to be tuned.The genetic algorithm cannot be applied since the computational demands are far too high.Therefore the tuning of the parameters has been performed on an empirical way.The calculation was performed with presumed parameters.The results (blood pressure waves) were examined and the parameters were modified in accordance.Blood pressure waves from Avolio [3] and Wang and Parker [18] were used for this process.

Results
After fine tuning the parameters of the viscoelastic branches a final calculation was carried out.The results of this simulation are treated as the healthy arterial blood flow.Blood pressure waves in the arterial branches A01 (ascending aorta), A37 (left external iliac) and A50 (left anterior tibial artery) are shown in Fig. 5.The curves represent one heart cycle.The quasi steady state was reached after a few (4-5) heart cycles, but the calculation was continued to eliminate also the minor transient effects.Therefore the analysis of the resulting curves was carried out 30 seconds after the initial state.Fig. 6 shows the blood pressure wave and the blood velocity curve in the A38 (right external iliac) branch.In this case two consecutive heart cycles were examined four cycles after the initial state.One can notice only slight differences between the two cycles, the quasi steady state is practically reached.The effect of arteriosclerosis was investigated as a next step.According to the advices of physicians the external iliac was examined in the first calculations.Therefore a stenosis was placed into the A38 arterial branch.The initial diameter of the stenosis was set to the initial diameter of the branch (8.0 mm).The diameter of the stenosis was then reduced stepwise until the value of 2.2 mm.The resulting changes in the blood flow were observed through these steps.A total of six calculations were carried out with consecutively reducing stenosis diameters.Fig. 7 shows the blood pressure wave and velocity curve in case of the least stenosis diameter.The volume flow rate as a function of the stenosis diameter is presented in Fig. 8.

Discussion
In the pressure curves of Fig. 5one can notice that the maximum of the amplitudes (i.e. the systole) is increasing as we move further away from the heart.At the same time the minimum (i.e.diastole) is practically unchanged.This phenomenon is in accordance with the experience of physicians.Also the average value of the blood pressure is increasing towards the peripheries.In the blood pressure curve of the aorta (A01) the so  called incisure point (dicrotic notch) is clearly visible just like on real measurement data.This also points at the agreeable accuracy of the arterial model.Blood pressure and velocity curves of the right external iliac (A38) are shown on Fig. 6.It should be noticed that during diastole the velocity is negative for a short period.Thus blood flow is reversed for a short time in this arterial branch.Fig. 7 presents the blood pressure and velocity curves for the same arterial branch (A38), however a stenosis with diameter of 2.2 mm is placed at the middle of the branch.The presented curves are observed right after the stenosis.The systolic pressure drops from about 122 mm Hg to 79 mm Hg.The change in the diastolic pressure is less: it drops from 77 mm Hg to 69 mm Hg.The peak of the velocity curve is also reduced: the maximum velocity is reduced from 0.52 m/s to 0.14 m/s.On the other hand the shape of the velocity curve is fundamentally different.No backflow occurs, in fact the velocity is positive during diastole.This is the reason why the volume flow rate is only reduced from 517 ml/min to 402 ml/min.A diameter reduction of 72,8 % only causes a 22,2 % drop in volume flow rate.This result was in accordance with the experience of physicians.

Conclusion
In the current paper a one dimensional arterial model network was created for modeling arterial blood flow.The viscoelastic effects of the arterial wall were taken into account with the so called Stuart model.Equations of the material model and the governing equations describing fluid flow were solved using the MOC.The so called Transient Simulator software was used to carry out calculations.The model network consists of 45 viscoelastic branches.The heart is modeled as a volume flow rate source, arterioles and capillaries are treated as linear resistances.Blood is assumed to be a homogenous, Newtonian liquid.The parameters of the Stuart model were tuned on an empirical way.Calculations were carried out on the resulting arterial model network.The results agree well with experiences of physicians.The effect of a stenosis in the right external iliac was investigated in the next step.The diameter of the stenosis was reduced stepwise, pressure waves and velocity curves were observed at the downstream end of the iliac.One could notice a significant drop in the systolic pressure.However there is a less reduction in the diastolic pressure.The general shape of the velocity curve changes remarkably.An interesting result is that reducing the diameter of the stenosis from 8.0 to 2.2 mm the volume flow rate of the blood reduces only 22 %.This result is confirmed by physicians.In the future the calculations should be compared to blood pressure measurements carried out on humans in order to validate the arterial model network.A mathematical method has to be developed for the fine tuning of the viscoelastic parameters.Once this method is functional, the arterial model can be adjusted to individual patients.The effect of arteriosclerosis can then be investigated in specific points of the model.The resulting pressure waves and blood velocity curves of the calculation are useful data in terms of the therapy.

Fig. 4 .
Fig. 4. Volume flow rate as a function of time -the boundary condition at the heart.

Fig. 5 .
Fig. 5. Blood pressure curves at the upstream endings of the arterial branches A01, A37 and A50

Fig. 6 .
Fig. 6.The effect of reducing the diameter on the volume flow rate (A38).

Fig. 7 .
Fig. 7. Blood pressure and velocity curves in the A38 branch with stenosis diameter set to 2.18 mm.

Fig. 8 .
Fig. 8.The effect of reducing the diameter on the volume flow rate (A38).