Mathematical Model of Arterial Hemodynamics , Description , Computer Implementation , Results Comparison

Analysis of basic parameters of a blood pressure in major arteries of human body is essential issue nowadays. It can help to select more appropriate way of surgical intervention for people who are ill with pathological changes of arteries, for example atherosclerosis. New mathematical model of blood pressure was developed to do such kind of analysis. The mathematical model is one-dimensional in terms of spatial coordinates. It contains verity of parameters that makes available to adjust this model to arterial bed of any exact case of any exact patient. Piece of software was implemented on the basis of the mathematical model. The software is aimed to build any particular part of arterial system and calculate parameters of blood pressure in it. Comparison of results gathered by developed software and by finite-element analysis software ADINA showed a slight difference. But the calculation time for the first one is a considerably less.


Introduction
The aim of the study was to address medical and social problems associated with the optimization of surgical treatment of blood circulation disorders in general and as a particular case of lower limb blood flow disorder.Statistics says that cardiovascular disease (CVD) is one of the leading causes of people disability and premature death in economically developed countries.Surgical reconstruction is widely used for the treatment of such diseases.However, there are still no objective indications which can help to choose a particular type of an implant or a shunt, its geometrical parameters, configuration and type of reconstruction.Because of this it is vital to develop a mathematical model, which should properly characterize the actual interaction of blood flow with the vascular wall.Such model should be fast (in terms of its implementation on a PC) and should have a set of parameters, which would allow adapting it to the particular patient.In addition, there is a necessity to develop simple user-friendly software for personal computer (PC), using the mathematical model.This software should allow quickly create a model of a part of a vascular system and calculate parameters of blood flow in this part.Such kind of software would allow simulating different variants of surgical intervention and would help to choose the most appropriate variant for every particular case and every individual patient in advance.

Mathematical Model
Existing mathematical models of blood flow do not have the necessary set of features and parameters.Because of this new mathematical model was introduced to match our requirements.The model is linear and one-dimensional.The main system of equations of the model can be solved analytically.In addition, this model takes into account the initial tension in the vessel and the angles between the arteries in a bifurcation or in a fork of other kind.

Main System of Equations for Blood Flow Dynamics Problem
Consider the system of equations that includes the following relations: simplified one-dimensional differential equation of a viscous incompressible fluid (Guljajev & Kossovich, 2001): there Q-volumetric velocity of a blood flow, P-pressure in blood, μ-blood viscosity, ρ-blood density,R-radius of a particular vessel, z-longitudinal coordinate, t-time; continuity equation which relates the volumetric flow rate Q with the radial displacement of the vessel walls w: (2.1.2) dynamic equations of motion of an axisymmetric circular cylindrical shell of the membrane theory (Pedli, 1983): where h-thickness of the vessel wall, u-longitudinal displacement of the vessel wall, K 1 and K 2 -coefficients of the compliance of tissues in the radial and axial directions, S and T-pulsation components of axial and transverse tensile force of the vessel, S 0 and T 0 -their initial values; relations of ideal elasticity of the vessel walls for general state of plane stress there E-elasticity modulus, v-Poisson's ratio; or well-known relations that take into account the anisotropy of the vessel wall (Lehnitskiy, 1977).
Thus, we obtain a closed system of six Equations (2.1.1)-(2.1.6)with respect to partial six unknown values u, w, Q, T, S, P. It is worth noting that there as in most papers blood is suggested as a liquid that has a constant viscosity and can be modeled as a Newtonian fluid (Berger & Jou, 2000;Tambasco & Steinman, 2003).At the same time Chen and Lu (2006) found that the difference in the values obtained for the case of Newtonian and non-Newtonian fluids do not exceed 10% for large vessels.In addition, Equation (2.1.3)-(2.1.6)allows taking into account compliance of the vessel walls, which is an important factor in modeling the flow of blood in the vessels (Rutten, 1998).Next step is to simplify somehow this system of equations.At first express pressure from (2.1.4)and we will have following formula Now it is possible to get rid of the pressure in relation (2.1.1)using this formula for pressure Substitution expressions for forces (2.1.5)and (2.1.6)in (2.1.8)leads to following . (2.1.9) After grouping terms with the equal derivatives (2.1.9)can be written as: The same transformation is applied to Equation (2.1.3)and it takes form of following formula (2.1.11) After grouping we have final equation (2.1.12) After all described transformations we have shorter system of equations instead of system of six equations for six unknown values: It is a closed system of three partial differential Equations (2.1.13)for three unknown functions . This system will be used as a basis for modeling the dynamics of blood flow in the arterial system.

Analytical Solution for the System of Equations
The unknown functions Q w u , , can be represented in the form of complex Fourier series: there T-the period of blood circulation.Next step is to substitute (2.2.1) in system of Equation (2.1.13)and then omit subscript k near unknown functions.After reduction by a factor t i k e  all equations of the system it can be written as Then we solve each equation of system (2.2.2) with respect to members containing higher derivatives of unknown functions.
In the second equation of the system we can get rid of the derivative 2 2 dz U d using the first equation.After that some simple mathematical transformations is applied and as a result following system of equations is obtained It is a system of three ordinary differential equations.It is worth noting that the unknown functions and coefficients of the equations are complex functions and values.Now we transform the system of three differential equations of second, third, and first-order (2.2.4) to a system of six first order differential equations.
For this purpose we introduce the auxiliary functions as follows: Substitution of these functions to the system of Equation (2.2.4) leads to the following result: Following constants are introduced for brevity: Assuming this constants system (2.2.6) can be written as Below is a matrix form of the system of Equations (2.2.8): -matrix of coefficients.
The system of differential Equations (2.2.8) or the matrix differential Equation (2.2.9) is the final relations that must be addressed.To solve the resulting system of equations we need to form characteristic equation (2.2.10) Calculation of determinant leads to a complex equation of sixth degree The characteristic equation is the sixth degree equation with complex coefficients.Since the equation contains only even powers of the variable λ it can be converted to an equation of third degree: Then every solution of (2.2.12) will correspond to two solutions of (2.2.11), . Thus, we get six eigenvalues of the system of Equations (2.2.9).We can find the eigenvectors for each of the six eigenvalues, by solving a following system of equations -eigenvector corresponding to the eigenvalue i  .Equation (2.2.13) can be written in a form of system of six algebraic equations An analytical solution of (2.2.14) gives the following expressions for the components of the eigenvectors The eigenvector for the eigenvalue λ i can be represented as follows: (2.2.16) As a result unknown functions are taken a form of following formulas: (2.2.17) Below are expressions for six unknown functions for system of Equations ( 2 The relations (2.2.18) are the solution of the original system of equations, but they are determined up to arbitrary constants of integration.For each artery it is needed to determine six arbitrary constants.

Boundary and Contact Conditions
Arbitrary constants appearing in relation (2.2.18) should be determined from the boundary and contact conditions.Under contact conditions are meant relations which are satisfied in a node where set of arteries are connected to each other (e.g.bifurcation or one to one connection).Such node in this article is called "contact node".Six conditions should be set for each vessel included in the part of vascular bed that is under consideration.The following set of relations can be taken as such conditions (assume that the part under consideration has only one beginning (input) artery and a number of ending (output) arteries): at the beginning of input artery should be set input volumetric blood flow velocity and fixation conditions: -at the output of ending arteries should be set relation between pressure and volumetric blood flow velocity and fixation conditions: there l-the length of the vessel.
-at contact nodes should be set following conditions: (total volumetric blood flow velocity of incoming and outgoing arteries is zero) , there n-number of arteries connected in the node, i l -circumference of cross-section of i-th vessel in the contact node.
Relations (2.3.1),(2.3.2) and (2.3.3)form a system of equations.There we have three conditions at the beginning and three conditions at the each ending of the part of arterial system that is under consideration.To make the system closed it is necessary to have three conditions for each of the arteries at each of contact node.Indeed, we have a balance equation for volumetric blood flow velocity, the two averaged equations of equality of longitudinal and transverse forces and n-1 equations for each pressure and both of displacements: 3 + 3(n-1) = 3n, i.e., for n arteries in contact node we have 3n equations, in other words 3 equations for each artery.This means that we have a closed system of algebraic equations to find all arbitrary constants for each artery included in the considered part of arterial system.To solve this system of algebraic equations first relation that describes input volumetric blood flow velocity should be decomposed into a complex Fourier series (2.3.1) .) ( Now we are trying to find solution for the pulsating flow sub-problem (only pulsating component of flow is under consideration) and because of this we can discard the constant term in (2.3.4).This term will be taken into account later during solving of steady-state flow sub-problem.Next step is to substitute into the boundary and contact conditions expressions (2.2.18).Then equate the terms at the same frequencies ω k and divide by After all these transformations we have following relations for the boundary conditions at the input: Conditions at the contact nodes should be transformed in the same way, but their form is depended on the configuration of particular contact node, the number of incoming and outgoing vessels, their mechanical and geometrical parameters.For a complete solution of the problem it is needed to solve steady-state flow problem.Eventually the final results should be the sum of correspondent parameters calculated for results pulsating flow sub-problem and for steady-state flow sub-problem.

Simplification of the Model (for Pulsating Flow Problem)
It is possible to simplify the solution of the problem by reducing the number of arbitrary constants of integration in each artery included in the part of the arterial system.This can be done for reasons of limiting in some respects.For example, if mass density of the material of the vascular wall tends to zero the two roots of equation (2.2.11) will also tend to zero.Such passage to the limit corresponds to neglecting the inertia forces of the vessel walls.In this case, we choose four eigenvalues at each artery, which does not tend to zero.Boundary and contact conditions can be taken like below: at the input: -at the outputs: -at contact nodes: . Thus, we have two conditions for each artery at the beginning and endings points, and two conditions for each artery at each contact node.In other words, we have 4 equations for each artery.That means that we have a closed system of algebraic equations for the arbitrary constants.
In the simplest case, only two eigenvalues of the four can be left.When longitudinal tensile force of the vessel walls S 0 tends to zero only two eigenvalues remain finite.Then the boundary and contact conditions can be taken as following: at the input: -at the outputs: at contact nodes: In this case we have two conditions for each artery to determine the two arbitrary constants which means that we have a closed system of equations.

Solution for Steady-State Flow Problem
This solution taking into consideration angles between arteries in contact nodes.At first we have to calculate the resistance to flow in contact nodes taking into account the angles between the incoming and outgoing arteries.As basic relations we take the dynamic contact conditions (Pedli, 1983).These equations were obtained from the conservation equations of a momentum of a continuum: We can leave in these equations only terms corresponding to the steady-state flow: Let's consider the case when there are only two arteries connected in the node: one is incoming and another is outgoing.For this we need just discard terms corresponding to the third artery.In addition should be taken in to account fact that because there are now only two arteries both of them have the same volumetric blood flow velocity ( ). Eventually we have following equations: Now we multiply first equation by cosα 2 and second one by sinα 2 and sum resulting equations.After such transformation we have next formula: . cos cos Next we consider electrodynamic analogy for outgoing artery (method of electrodynamic analogy for blood flow analysis was first used in study (Guljajev & Kossovich, 2001) for the case of internal mammary artery).Electrical schema for this analogy is shown on the Figure 2.5.1.There P 11 -pressure at the end of the incoming artery, P 20 -pressure at the beginning of outgoing artery, P 21 -pressure at the end of outgoing artery, R y -resistance of contact node, R π -Poiseuille resistance of the outgoing artery, Q 2 -volumetric velocity of blood flow.Pressure is an equivalent for electric potential, resistance is an equivalent for electric resistance and volumetric velocity is an equivalent for electrical current.Due to the fact that this scheme represents sequencing circuit current in each point of it is the same.According to the Ohm's law we have following relations: In terms of Figure 2.5.1 Formula (2.5.1) for part with resistance У R can be written as: (2.5.3)Substitution of (2.5.3) into (2.5.2) leads to: Then we transform Formula (2.5.4).Below is a result of this transformation.After substitution (2.5.6) into (2.5.5) we express R y from resulting relation due to the small pressure drop in the blood vessels relative to normal atmospheric pressure.
To calculate parameters of steady-state flow we need to combine a system of equation for whole part of arterial system that is under consideration.This system should contain following equations: for each artery relation of volumetric velocity of blood flow and pressure (analogy for relation of electrical current and potential) according to Formula (2.5.2) (2.5.8) there k-number of incoming artery in contact node; i-number of outgoing artery; 1 , k P -pressure at the end of k-th artery, in other words, pressure right before contact node; 1 , i P -pressure at the end of i-th artery.For the first (input) artery of considered part of arterial system we have relation .5.9) there k i -set of numbers of arteries connected in the node -at the beginning of input artery we set initial volumetric velocity of flow Q 1 = Q 00 , there Q 00 -constant term of the series (2.3.4).
-at the endings of output arteries we set pressure P jk,1 = 0, there j k -set of numbers of ending (output) arteries for the part of arterial system that is under consideration (there P i,0 and P i,1 -pressure at the beginning and at the ending of i-th artery correspondingly).
Thus we have closed system of equations to calculate parameters of steady-state flow in the arterial system taking into account the angles between the vessels in contact node.

Example of System of Equations for Steady-State Flow Problem
There is shown on a particular example how to combine such system of equation.Example arterial system part for the example is shown on a Figure 2.6.1.The part of arterial system under consideration consists of three arteries and one contact node.The contact node has one incoming artery and two outgoing arteries.For such vessels bed we have electrodynamic analogy that is shown on the Figure 2.6.2.Next step is to analyze this schema according to Formula (2.5.8).In first part we have pressure drop (potential) equal to P 1,0 -P 1.1 .Resistance to flow in this part is equal to Poiseuille resistance due to the fact that the first part is an input and does not come out from any contact node: For second part drop of pressure is equal to 1 , 2 1 , 1

P P 
. Resistance in this case consists of two components, Poiseuille resistance and resistance of contact node and equal to . Relation for volumetric velocity can be written as: (2.6.2) Similarly we build relation between volumetric velocity and pressure for third part: (2.6.3) Resistances of contact node iУ R are defined by (2.5.7).Poiseuille resistance can be calculated using following formula (Guljajev & Kossovich, 2001): there μ-viscosity coefficient of blood, l i -length of vessel, r i -vessel radius.
According to (2.5.9) we have following relation for contact node: (2.6.4) We need to set input volumetric velocity: (2.6.5)Also we need to set pressure at endings: (2.6.6) Finally we have a system of algebraic Equations (2.6.1)-(2.6.6) that can be written as: This is a closed system of seven algebraic equations with seven unknown values Q 1 , Q 2 , Q 3 , P 10 , P 11 , P 21 , P 31 .This system can be solved using different mathematical approaches e.g. using Gauss method.We can find all needed parameters of steady-state flow by solving this system of equations.Thus, the problem of the pulsation of blood flow is finally resolved.The final results are the sum of the steady-state and fluctuation components.

Determination of Excessive Blood Volume in the Vessels
One of the most important parameter of the blood circulation is an excessive volume of blood in vessels.To determine this parameter we consider equation of continuity for i-th vessel. There After simple transformation we have following formula: Integrating by dz the last relation, we obtain expression Integrating this by equation by time give us formula for excessive volume of blood in vessel in exact period of time (2.7.2)

The Software Program
New piece of software has been developed to perform the actual calculations using a mathematical model described above.The software program is designed to simulate blood flow in the considered part of the human arterial system.The software program allows you to graphically build vessels bed to set the parameters of blood and each vessel separately.Algorithm assemble at run-time systems of equations for given vessels bed.After all systems of equations are assembled algorithm solve them and output solution as plots of different kind into multi-windows user interface.Below are lists of input and output parameters of the software program.
Input parameters:  Geometry of arterial system (vessels bed configuration).


Mechanical parameters of blood (density, viscosity). Mechanical characteristics of vessels (Young's modulus, Poisson's ratio, the initial tension). Volumetric velocity of blood flow at the beginning of input artery -Q (t).


The parameters relating the volumetric velocity of blood flow with a pressure at endings Output data:  Blood flow pressure in arterial bed.


Volumetric blood flow velocity in arterial bed. Linear velocity of blood flow (average over the cross section).


Excessive volume of blood in arteries.
Results are represented as plots of relationship between different parameters of blood flow and longitudinal coordinate and time.There is ability to get two-dimensional and three-dimensional plots.Two-dimensional plots show the change of parameter along one of variable with the second one fixed.Also it is possible to view animated graphics where each frame shows plot for particular value of fixed variable.One more feature is that it is possible to view plot for particular artery or a set of plots for every artery in system in one window.Screenshots of the software program user interfaces are shown on Figures 3.2.1-3.2.5.

Comparison of the Results
To understand whether a particular mathematical method or program gives adequate results or not it is useful to compare its results with experimental data or with data obtained by any other methods or other programs.In this research we compare results obtained using described model and software program with a results obtained using the finite element method.ADINA System 8.2 software program was chosen as a finite element method implementation.Three-dimensional model was built in a CAD system Solid Works 2007.

Input Parameters (Characteristics of the Model)
Bed of the femoral artery of the author was chosen for calculations (as object for simulation).All parameters of the arterial bed were obtained by duplex ultrasound device Toshiba Xario.A geometric model of the arterial system of the femoral artery has the parameters shown in  (Begun & Shukeilo, 2005), Poisson's ratio v = 0.5 (Vessel wall is incompressible). The boundary conditions.The inputs and outputs cross sections are rigidly fixed.The inner surface is the surface of contact with blood.No-slip conditions are selected on the contact surface (the fluid velocity matches the velocity of the wall on this surface).
 Finite element is a tetrahedron with edge length 0.001 m.As a result after partitioning we have about 74 791 finite elements.

Blood model has following mechanical parameters:
 Density: (viscosity can be measured using rotational viscometer at velocity equal to 100 c -1 ). The boundary conditions.The lateral surface is the surface of contact with the vessels wall.No-slip conditions are selected on the contact surface.Front surfaces at outputs are considered as free surfaces.This means that pressure here is equal to zero.At input front surface we set relationship between initial blood flow velocity and time (

Calculated Results
As a result of calculations we have plots of dependence of average velocity on time at input and output cross-sections of the arterial system that is taken into consideration.This plots are shown on the Figures 4.2.1-4.2.4.Next step is to perform the same simulation using the mathematical model and the software program described in this article.For this purpose we use the same input parameters (same part of arterial system, same blood and vessels parameters) as for simulating using finite element method.Results of calculation can be found on the Figures 4.2.5-4.2.8

Compare Results
To make comparison of results more illustrative results are exported into Microsoft Office Excel.Then we build diagrams for the results acquired for both methods in the same cross section on the same coordinate plane.Result of such data massaging are shown on Figures 4.3.1-4.3.4.These diagrams show that the results obtained by finite element method and method described in this article are close to each other.The maximum deviation between results of these two methods we have at the moment of systole and this difference does not exceed 10%.Elsewhere in the diagrams curves differ slightly or even coincide.The curves in Figure 4.4.1 coincide completely due to the fact that it represents input cross section and we set up flow velocity here manually as an input parameter for both of methods (it is a boundary condition).Thus, the developed mathematical model and software system can be used to analyze the overall, global state of the periodic laminar viscous flow in elastic tubes and as a variant of the blood flow in the arterial bed.However, since that is one-dimensional mathematical model, it does not allow us to analyze the distribution of a parameter in a small neighborhood of a point (or aria) with a strong geometric or physical heterogeneity, for example in the small neighborhood of the fork.Also this model allows calculating the flow velocity averaged over the cross section, but it does not allow obtaining the velocity profile in a particular section.Such problems are easier to solve using the finite-element simulating.2) It is shown that results obtained using one-dimensional model differ slightly from the results of calculations obtained using three-dimensional model.However, since that is one-dimensional mathematical model, it does not allow us to analyze the distribution of a parameter in a small neighborhood of a point (or aria) with a strong geometric and physical heterogeneity, for example in the neighborhood of the vessels fork or near atherosclerotic plaques.Also, the model does not allow analyzing the velocity distribution over the cross section and does not account bending of the vessels.
3) The new software program was developed.This software program implements described mathematical model.This software program can simulate a wide range of the vascular system configurations and easily customized to a specific case.
4) The described mathematical model and software program can be regarded as a basis for further clinical studies to substantiate the selection of the method and option of reconstruction including type and shape of the plastic material to suit the individual characteristics of arteries of each patient.
in formulas here insignificantly differ from the normal atmospheric pressure we can

Figure 3
Figure 3.2.1 User interface intended to build arterial bed and set arteries parameters Figure 4.1.8)(m/c -c). Finite element is a tetrahedron with edge length 0.001 m.As a result after partitioning we have about 73 895 finite elements. Initial volumetric blood flow velocity at input is selected to be as it is shown on Figure 4.1.1.

Figure 4
Figure 4.1.1.Plot of relationship between velocity and time (m/c -c).Double period of pulsation

Figure 4
Figure 4.2.1.Plot of dependence of velocity on time in the input cross-section of the common femoral artery

Figure 4
Figure 4.2.4.Plot of dependence of velocity on time in the output cross-section of the peroneal trunk

Figure 4
Figure 4.2.7.Plot of dependence of velocity on time in the output cross-section of the popliteal artery

Figure 4
Figure 4.3.1.Plot of dependence of velocity on time in the input cross-section of the common femoral artery

Figure 4
Figure 4.3.2.Plot of dependence of velocity on time in the output cross-section of the deep femoral artery

Figure 4
Figure 4.3.3.Plot of dependence of velocity on time in the output cross-section of the popliteal artery