Размещено на http://www.allbest.ru/

Government of Russian Federation

National Research University

Higher School of Economics

Faculty of Computer Science

Applied Mathematics and Information Science Program

Diploma work

«Calibration of 1-D blood flow model for real patient»

Student:

Group 122

E.A. Sakharova

Research advisor:

T.K. Dobroserdova

Moscow 2016



Table contents

1. Introduction

2. Model

3. Test with silicone model

4. Results of tests with silicone graph

5. Calculations on real patient vascular tree

6. Results of real-patient graph calculations

7. Conclusion

8. Preferences



1. Introduction

Worldwide major causes of death are cardiovascular diseases such as atherosclerosis, thromboembolism, aneurysm and many others. An important problem of modern medicine is to create effective methods of diagnostics, treatment and prevention of abovementioned diseases. Mathematical modeling and numerical calculations of blood flow in networks of human vessels play a significant role in medicine. Without a doubt, a profound understanding of pressure and flow pulse wave propagation in the cardiovascular system, the impact of disease and anatomical variations on these propagation patterns can provide sensible information for medical treatment.

Already in the XVI century, scientists achieved the great results in research of heart and blood flow. British medic William Harvey is considered as a founder of this study. He substantiated the difference between breath and pulse, described the heart work process and physiology of vessels. In addition, the study of blood circulation might not exist without Bernoulli and Poiseuille works. Flow study unites both physicists and medics.

However, modern mathematical modeling of human cardiovascular system have undergone a vast development over the last few decades. Researchers have theoretical insight of cardiovascular system. Varieties of clinics help with necessary measurement and data. Thanks to rapid grows of computational technologies quick and accurate calculations of big data are no longer a problem. Specifically, various deterministic models may give the presentation of the full arterial tree or target specific sites.

Also there is abundant literature concerning the problem of modeling of arterial system. The present researches involve numerical calculations using governing equations from zero-dimensional to three-dimensional. In more details, 0-D model uses in electrical circuits [14], 1-D model is the most popular and convenient, because the problem statement is simple and calculations are not very time-consuming. 2-D and 3-D models usually applied for blood simulation in a small part of vascular network (one or several vessels). It is necessary to obtain additional unknown parameters in assuming wall elasticity. example, the article [15] provides a results and discussions about use of 1-D and 3-D models. Moreover, they found that 1-D model is a reasonable representation of 3-D system and can be used to estimation of various conditions. Calculations on 1-D model less complicated and their results can be applied in 3-D formulation.

Majority of articles base their calculations on the various 1D mathematical models. These differences involve boundary conditions on the end of each vessel in assumptions on the vessels' wall material. So called equations of state help to take into account such properties. According to [19] the detailed description of vessels' structure is given. In the same publication authors provide different dependance between pressure and cross-sectional area using empirical curves. However, the work [6] studied the role of visco-elasticity in explaining the divergences of the results. The small reductions in the comparative errors of pressure and flow rate suggested that the effect of visco-elasticity on pressures and flow was secondary.

As arterial as venous system can have different graph representation. Of course, it is impossible to measure all vessels (including capillaries) and construct a graph of full blood system. It is clear that the bigger arteries the model involves the better. Consequently, calculations on the big graph become more complicated. That is why, at first, there were a lot of calculations on small graphs consisting about 30-55 vessels. Now there are a lot of calculations on the bigger graphs. The paper [15] provided calculation on large and small arteries separately and gave more detail information about the form of the waves change and geometry. In priority, they proved that the structured tree model can predict flow and pressure not only in the large but also in the small arteries. Some scientists adjust only the head graphs [2, 6]. Except this, in the work [7] graph of 2142 arteries was constructed. Authors configure all parameters in the hemodynamics model. It has patient specific description and predictive capabilities.

A modern problems in hydrodynamics faces with the problem of solution systems of partial differential equations. The book [17] contains derivation and a detailed description of grid-characteristic methods. As mentioned, the numerical methods of solving the system of differential equations are very time-consuming. Different researchers use various schemes in attempt to choose suitable one in the face of this tricky task. The article [1] may be helpful in this kind of job. Scientists from the United Kingdom and Brazil reviewed the six popular numerical methods and compared them with each other. As the result, all six schemes considered in this study have been assessed in a series of controlling tests with an increasing degree of complexity, for which theoretical, numerical or artificially created pulse waveforms are available. The research confirmed a good agreement among all numerical schemes and their ability to solve the equations.

Finally, all these studies provide additional support for the use of 1D reduced-order modeling to accurately simulate arterial pulse wave with a reasonable computational cost. The majority of studies use graphs of an abstract patients and do not take into account individual patient features. This work is applied to research and customize parameters in the model using real patient data. Moreover, the choice of equation of state, boundary conditions and numerical schemes is still open. The goal of this study is to apply hydrodynamics model to real-patient graph and discuss results with doctors.

2. Model

The human organism has systemic and pulmonary circulation (Fig.2.1). The systemic circulation serves to deliver to all organs and tissues of the body needed amount of nutrients and oxygen, while the pulmonary circulation serves to enrich the blood with oxygen in the lungs.



Fig.2.1. Systemic and pulmonary circulation.

Circulatory system is divided into two connected parts: venous and arterial (Fig. 2.2).

Fig. 2.2. Venous (left) and arterial (right) circulatory systems.

Fig. 2.3. Example of graph representation of 61 main arteries of human organism [1].

In the model of blood circulation vascular system can be represented as a graph. In particular, vessels corresponds to edges of this graph, nodes - docking points (for example Fig.2.3). Obviously, vessels are considered as elastic tubes, each of them has specific properties. These tubes can be represented by a surface of revolution (Fig.2.4) in which ratio of the average radius to the length is small.

Blood is supposed as a viscous incompressible Newtonian fluid and its flow in each one-dimensional region can be described by the laws of conservation of mass (1) and momentum (2) :

, (1)

(2)

where

velocity averaged over cross-sectional area;

S - cross-sectional area of the vessel;

- transmural pressure (difference between internal and external pressure);

t - time, , T - calculation time;

- coordinate along the vessel, - length of the vessel;

- blood density;

- specified functions.

In detail, function may define the blood inflow/outflow, and function - the impact of external forces (friction force, gravitation) [6, 15, 18, 19].

Fig. 2.4. Example of vessel representation, where r - radius, L - length.

The system (1)(2) closes by equation of state, that characterizes elastic properties of the walls of vessel [5, 9, 15, 18, 19]. In scientific literature it is known various ways to its' formulation.

For example,

1. The derivation from the equations of equilibrium under the assumption of linear elasticity for the wall material:

[2, 5, 6], (3)

where

, - average wall thickness, E- Young's modulus;

- the Poisson's ratio, - cross-sectional area when ;

2. Empirical function

[9, 18, 19], (4)

where - speed of pulse wave propagation. In this case equation of state takes into account the non-linear wall material.

More complicated equations of state were derived in the assumption of nonlinear vessel's wall material properties [5].

The system (1)(2) has hyperbolic type, consequently, characteristic curves impose compatibility conditions on the ends of each vessel. For the correct formulation of the problem, it is required one additional conditions on the end of each vessel. So at the junction of the branches it is necessary to define the system of equations consisting laws (5)(6):

the conservation of mass law: (5) and

the law of Poiseuille flow:

(t, ) (6),

where - input/output flow,, n - the number of vessels, - pressure in vessel №i, - resistance, imitated microcirculation №i, --pressure in node;

It is common to use Poiseuille law as a boundary condition. However, it has the unknown variables which are very complicated to define. That is why we use alternative forms of boundary conditions instead (6) using (5). They can be represented in different forms a) - c):

(continuity of total pressure);

(continuity of pressure);

c) for bifurcations: [16],

where index «1» corresponds to mother branch, and «2» и «3» - to daughter branches;

, i = 2,3; the angle of inclination of branch №i relatively the parent branch.

For system closing, this task required conditions on terminal nodes, which can be represented in following ways:

i) (Poiseuille condition)

, n - the number of vessels, - pressure in vessel №i, - resistance, imitated microcirculation №i, -- terminal pressure;

ii) Imitation of electrical circuits [14];

iii) Impedance boundary conditions [3, 12].

This conditions also came from the theory of electrical circuits. Obviously, that it is desirable to take into account microcirculation on the ends of the main arteries. Consideration of these vessels in the structure of the graph makes calculations very complicated.

Consider, that the value is constant for every vessel in the graph. In addition, every terminal vessel has the following structure corresponds to microcirculation (Fig.2.5): each vessel bifurcates into two child vessels unless one of their radius is bigger than . Radiuses of daughter vessels calculates as . Assume, that bifurcations have the same structure as on Fig.2.6.

In the model of blood circulation impedance is used as a replacement of resistance. The simplest case includes use of Darcy's law P=RQ, where R - resistance value. It is possiable to apply Laplace transform .

System (1)(2) can be rewritten in terms of . This system has solution evaluated at x=0. This solution defines through defines through It means, that impedance on the beginning of the vessel defines through impedance on the end of the vessel. On the each smallest vessel (where ) we define terminal impedance and can calculate for this vessel (Fig. 2.7). Impedance calculates recursively:

Step 1: Impedance on the terminal branches at the end of vessel:

Step 2: Impedance on the terminal branches at the beginning of the vessel: ;

Step 3: Impedance in bifurcations. It corresponds to the impedance at the end of parent vessel.

In the bifurcations impose the condition: ,

where - parent vessel impedance, - daughter vessel impedance.

Reffering to [3] we obtain relation for each bifurcation and impedance satisfies. It means that now equation (7) can be transformed corresponding time domain in form in the assumption of periodic solution:

.

Next, recieve descrete solution for pressure at n-th time step using inverse Laplace transform:

,

where - impedance weights and flow at the j-th time step, - terminal pressure.

iv) Moreover, the closure may be obtained in different way: use of venous and arterial trees. There is a drawback, that it is necessary to define docking conditions of these two trees. For example, in this case we can use condition i) but the problem of tuning resistances it still actual. Basically, this approach has a complicated problem of additional calibration of venous tree.

In this study a numerical model Simakov-Kholodov, MIPT was used for the calculations [9]. There are a lot of numerical methods for solution [17]. Grid-characteristic numerical method used here is a good alternative for Lax-Wendroff method. The first order scheme is conservative, monotone and contains essential dissipative terms for definition of non-smooth solutions. No shock waves appeared in the system, the coefficient of hybridity was selected in such a way that the scheme could be of the second order.

Fig. 2.5. Example of vessel structure on the terminal branches.

Fig. 2.6. Vessel bifurcation.

Fig. 2.7. Impedance algorithm.

3. Test with silicone model

Depending on the using boundary conditions and equation of state, the model of blood flow can have a number of unknown parameters: the wall thickness, the angles in the bifurcations, resistances, speed of pulse wave propagation. These parameters are specific to each patient, it is difficult to measure them, so we need to understand is it possible to avoid their use without loss of quality of the numerical solution.

Referring to the article [5] authors provided an experiment with simulation of fluid flow in human arterial network consisting of larger vessels. In particular, the physical model consisted of 37 interconnected silicone tubes. The tubes were attached to the pulsatile pump, which was providing a periodic input flow (Fig.3.1). In addition, the special fluid simulated the blood run on tubes. Sensors, installed on several tubes, measured tested parameters: pressure and flow. Tubes' parameters and material properties are known from the article. Authors also presented numerical calculations on the graph (Fig.3.2) corresponding to the model of silicone tubes. The calculation results have a small deviation from the experimental measurements.

Fig. 3.1. Experimental model of main human arteries.

In this work, we repeated the calculations on the same graph (Fig.3.2) using the parameters and conditions described in [5]. The article provides numerical and experimental curves of 9 “control” vessels (defined on Fig.3.2).

At the inlet of the ascending aorta, we take the following flow rate (Fig. 3.3). This inflow quite accurate matches with real inlet waveform.



Fig. 3.2. Graph, that corresponds to experimental model

Fig. 3.3. The experimental flow rate at the ascending aorta (vessel №1) taken from [19] t - time of one flow rate period, Q - flow, period T = 0.827.

Boundary conditions on terminal vessels defined as reservoir overflow of liquid. This reservoir has a constant volume, like human body.

Therefore, , where - the flow rate at the outlet of the 1D terminal branch, - the pressure at the outlet, - hydrostatic pressure- known resistance at the end of the silicone tube.

Known values:

, , ;

, - friction force, viscosity coefficient;

, c -speed of pulse wave propagation;

E = 1.2 MPa - the Young's modulus, .

To assess the significance of the model parameters and to improve similarity the numerical solution and experimental measures, the problem was calculated with different equations of state and the boundary conditions at the nodes.



4. Results of tests with silicone graph

Experiment 1:

We tested 1st and 2nd order schemes realized in Simakov-Kholodov model and compared pressure (Fig. 4.1) and flow (Fig. 4.2).

Fig. 4.1. The flow comparison for different numerical schemes. 1 - experimental curve; 2 - 1st order; 3 - 2nd order; 4 - numerical curve from article [5].

Both schemes gives the similar results and corresponds to experimental data. First order scheme is monotonous and smooths the solutions. Second order scheme is not monotonous that is why solutions have more fluctuations. Basically, second order is more time-consuming and we decided to apply first order for the next calculations.



Fig. 4.2. The pressure comparison for different numerical schemes. 1 - experimental curve; 2 - 1st order; 3 - 2nd order; 4 - numerical curve from article [5].

Experiment 2:

According to the admission that averaged cross-sectional area is constant on each vessel, obtained numerical pressure (P) and flow (Q) using boundary conditions item a) - c). The typical comparisons of results is shown on Fig.4.3 and Fig.4.4.

OX axis corresponds to time of one cycle (0.827 sec);

OY axis corresponds to flow in ml/sec (Pic.4.3) or pressure in kPa (Pic.4.4).

It is clear from the Fig.4.3 and Fig.4.4 that the results for boundary conditions are nearly the same. However, the profile of pressure and flow exactly reproduces experimental curves. We conclude that there is no need to use bifurcation angles. That is why we decided to provide calculations with one type of boundary conditions a) next.



Fig. 4.3. The flow comparison for different types of boundary conditions. 1 - experimental curve; 2-4 - a)-c) conditions; 5 - numerical curve from article [5].

Fig. 4.4. The pressure comparison for different types of boundary conditions. 1 - experimental curve; 2-4 - a)-c) conditions; 5 - numerical curve from article [5].

Experiment 3:

This experiment compares the results in case 1 and case 2 with experimental results (Fig. 4.5 and Fig. 4.6).

case 1: equation (3) with constant cross-sectional area

case 2: equation (4) with constant cross-sectional area.

It is clear from the plots, that both equations of state give the results closed to the experimental measures. For real applications the type of equation of state can be chosen in condition of known parameters.



Fig. 4.5. The flow comparison for different equations of state. 1 - experimental curve; 2 - eq.(3); 3 - eq.(4); 4 - numerical curve from article [5].

Fig. 4.6. The pressure comparison for different equations of state. 1 - experimental curve; 2 - eq.(3); 3 - eq.(4); 4 - numerical curve from article [5].

Experiment 4:

After results received from Experiments 1 and 2, we made an attempt to take into account variable cross-sectional area. In this case, calculations used different computational grid steps: 2 cm, 1 cm, 0.5 cm and 0.25 cm. It was established, that during the grinding less than 1sm the numerical solution almost does not change and corresponds to the experimental results. A more coarse mesh gives a considerable deviation from the “reference” values (to which the numerical solution converges with the refinement of the grid). According to these results, Fig.4.7 and Fig.4.8 shows results for variable cross-sectional area with optimal grid steps.

These results demonstrated that for applications where we need to reproduce pressure flow profiles it is possible to use constant cross-sectional area instead of various without accuracy solution losses. At the same time variable cross-section can be significant if we need to calculate reflected waves.

Fig. 4.7. The flow comparison for different cross-sectional area. 1 - experimental curve; 2 - variable cross-sectional area; 3 - constant cross-sectional area; 4 - numerical curve from article [5].

Fig. 4.8. The pressure comparison for cross-sectional area. 1 - experimental curve; 2 - variable cross-sectional area; 3 - constant cross-sectional area; 4 - numerical curve from article [5]. 

5. Calculations on real patient vascular tree

Segmentation

With the segmentation of images of computed tomography (CT) of the individual patient obtained three-dimensional region of the vascular bed. This process includes filtering the input image: bones, air and tissues separation [11]. Then extraction of vessel surface is produced. By extracting the centerlines of the three-dimensional vascular network there were reduced to one-dimensional graph. In result 1-D vascular tree described by a set of edges and nodes.

Parameters

The geometry and location of the vessel are individual. For every vessel we know length, cross sectional area, and for bifurcation - the angles. This graph can consist not only of the major arteries, but also take into account the small vessels.

We have a graph of arterial network of a particular patient (Fig.5.1 painted in Paraview program). This graph consists of 71 head branches (Fig.5.3-5.4) and connected 192 body branches (Fig.5.2). Lengths and radiuses of all vessels are known.

Unfortunately, it is impossible to measure all necessary parameters for every patient. Nevertheless, that some parameters can be defined depending personal features of human (age, chronic diseases, bad habits). Some of them such as Young's modulus or elasticity also can be found in medical atlases according mentioned features. However, resistances and wall thicknesses are still not found.

Model

Previously, we studied different boundary conditions, equations of state and convergence of two numerical schemes. In this part of work we apply these results and use different methods for graph calibration.

Mathematical model in this task is following:

1) system of (1)(2) equations with (4) equation of state (elasticity parameters for all vessels were found in medical literature);

2) the continuity of total pressure boundary conditions on internal nodes (no additional parameters are needed);

3) inflow defines as empirical curve on Fig. 5.5;

4) graph consists of 263 arteries was reconstructed for real patient Fig. 5.1.

Calculations provided using numerical model Simakov-Kholodov with 1st order numerical scheme (Ch. 2).

Fig. 5.1. 3D real patient graph representation of all 263 vessels.

Fig. 5.2. 2D graph structure of body arteries. 

Fig. 5.3. 3D representation of head arteries.

Fig. 5.4. 2D graph of head arteries.

Fig. 5.5. The empirical inlet flow Q.



6. Results of real-patient graph calculations

This part of study describes test with different boundary conditions on terminal branches.

Experiment 1: Nonreflecting boundary conditions [8] (condition on the second characteristic curve).

It was developed by Hedstorm, fully explained in [13] and has lot of applications in researches.

The main advantage of these conditions that there is no need to find additional parameters in hydrodynamic model. However, there is no blood mass control. In our experiment it leads the effect of ever-increasing pressure (Fig. 5.6). That is why these conditions are required to supplement.

Fig. 5.6. Flow and pressure in vessel 1.

Experiment 2:

On the terminal branches we regulated the inflow as follows:

Define Qin - input flow from heart (Fig. 5.5). We specify flow on terminal branches proportionally the diameter of terminal vessel.

, where m - the number of terminal vessels, - diameter of i vessel.

* Qin , i = 1..m.

Fig. 5.7 shows the results of flow in the basic arteries in head.

Fig. 5.8 provides results of flow in the two body arteries (№86, 117) and four atreries in legs (№ 160, 176 - right leg, № 236, 248 - left leg).

Fig. 5.7. Results of flow in control vessels in the head.

Some results are similar to medical measurements. However, the accuracy is still in question. Future work requires to use of another terminal conditions.



Fig. 5.8. Results of flow in control vessels in body and legs.



Conclusion

To summarize, this paper gives the basic concepts of human cardiovascular system, defines the system of partial differential equations to determine the process of blood circulation. Numerical model was based on the grid-characteristic method. Therefore, the developing method was applied for calculations of general variables.

In the first part of work we modelled the physical experiment with fluid flow in the silicone network. The blood flow model was verified with this benchmark. In addition, we tested various types of node boundary conditions, equations of state, variable and invariable vessel cross-sectional area, numerical methods of 1st and 2nd order to understand, which modifications of blood flow model are appropriate for patient-specific simulations. Consequently, for modeling of real-patient hemodynamics the following conditions were chosen:

· the continuity of total pressure boundary condition, because it does not involves additional parameters;

· equation of state (4), where elasticity can be found in medical atlases;

· invariable cross-sectional area, because calculations are significantly quicker than variable cross-sectional area, but result pressure and flow waveforms are very closed in both cases.

The second part of work includes calibration of real patient vascular tree. The graph was obtained from patient CT imaging. The different conditions on terminal branches were tested. All results are discussing with doctors.



Preferences

1. Alastruey J., Boileau E., Donders W., Eikeland Fossan F., Hellevik L., Huberts W., Mьller P., Nithiarasu P., Willemet M. A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling. In: Int. J. Numer Method Biomed Eng. 2015 Oct;31(10). doi: 10.1002/cnm.2732.

2. Alastruey J., Byrd S., Parker K., Peiro J. Modelling the circle of Willis to assess the effects of anatomical variations and occlusions on cerebral flows. In: Epub 2006, 40(8): pp.1794--805.

3. Alastruey J., Chowienczyk J., Epstein S, Willemet M. Reducing the number of parameters in 1D arterial blood flow modeling: less is more for patient-specific simulations. In: Am J Physiol Heart Circ Physiol 309: H222-H234, 2015.

4. Alastruey J., Figueroa C.A., Xiao N. A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models. In: Int. J. Numer. Meth. Biomed. Engng. 2014, pp. 204--231.

5. Alastruey J., Joaquim P., Koen M. Pulse wave propagation in a model human arterial network: Assessment of 1-D visco-elastic simulations against in-vitro measurements. In: J. Biomech, 40(15), 2007, pp. 3478--86.

6. Alastruey J., Khir A., Matthys K., Parker K. Segers P., Sherwin S., Verdonck P. Pulse wave propagation in a model human arterial network: Assessment of 1-D visco-elastic simulations against in vitro measurements. In: J. Biomech. 2011 Aug 11; 44(12): pp. 2250--2258.

7. Aurelio M., Blanco P., Lemos P., Passos R., Watanabe S. An Anatomically Detailed Arterial Network Model for One-Dimensional Computational Hemodynamics. In: IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 2, 2015.

8. Bo-Wen Lin, Pong-Jeu Lu High-resolution Roe's scheme and characteristic boundary conditions for solving complex wave reflection phenomena in a tree-like arterial structure. In: Journal of Computational Physics 260, 2014, pp. 143-162.

9. Dobroserdova T., Salamatova V., Simakov S., Vassilevski Yu. Blood flow simulation in atherosclerotic vascular network using fiber-spring representation of diseased wall. In: Math. Model. Nat. Phenom. Vol. 6, No.5, 2011, pp. 333--349.

10. Fukasaku K., Liang F., Liu H., Takagi S. A computational model study of the influence of the anatomy of the circle of willis on cerebral hyperperfusion following carotid artery surgery. In: Liang et al. BioMedical Engineering OnLine 2011.

11. Gamilov T., Ivanov Yu., Kopylov P., Simakov S., Vassilevski Yu. Patient Specific Haemodynamic Modeling after Occlusion Treatment in Leg. In: Math. Model. Nat. Phenom. Vol. 9, No. 6, 2014, pp. 85--97.

12. Guerrisi M., Karaja T., Kokalari I. Review on lumped parameter method for modeling the blood flow in systemic arteries. In: J. Biomedical Science and Engineering, 2013, pp. 92--99.

13. Hedstrom G. Nonreflecting boundary conditions for nonlinear hyperbolic systems. In: Journal of Computational Physics 30, 1979, pp. 222-237.

14. Himeno R., Liang F., Liu H., Takagi S. Multi-scale modeling of the human cardiovascular system with applications to aortic valvular and arterial stenoses. In: Med Biol Eng Comput (2009) 47:743-755 DOI 10.1007/s11517-009-0449-9.

15. Kim W., larsen J., Nadim A., Olufsen M., Pedersen E., Peskin C. Numerical Simulation and Experimental Validation of Blood Flow in Arteries with Structured-Tree Outflow Conditions. In: Annals of Biomedical Engineering, Nov. 2000, Vol. 28, Issue 11, pp. 1281--1299.

16. Lamponi D. One dimentional and multiscale models for blood flow circulation. PhD thesis, EPFL, Lausanne, Thesis n.3006, 2004.

17. Magomedov K., Kholodov A. Grid-characteristic numerical methods. In: Nauka, Moscow, 1988 (in Russian).

18. Mukhin S. Mathematical modeling of hydrodynamics. In: Moscow, 2008 (in Russian).

Размещено на Allbest.ru