modelling of the circulation

(1)

Minimal haemodynamic

modelling of the circulation

P.C.I. Spelde

Master Thesis in Applied Mathematics

April 2008

(2)

(3)

Minimal haemodynamic

modelling of the circulation

P.C.I. Spelde

First supervisors: A.E.P. Veldman and G. Rozema Second supervisor: A.J. van der Schaft

External supervisor: N.M. Naurits

Institute of Mathematics and Computing Science P.O. Box 407

9700 AK Groningen

The Netherlands

(4)

(5)

Abstract

The knowledge of the flowstructures in the human arteries is limited. The medical staff have the wish to have a better side to this phenomenon. In a specific mathematical research of the flow through the carotid bifurcation there is attention for this problem. To make it possible to do this research a mathematical model of the whole cardiovascular system (CVS) is needed.

Models found in literature simulate specific areas of the CVS while others are either overly complex, difficult to solve, and/or unstable. This thesis develops a minimal model with the primary goal of having the possibility to reflect accurately a small part of the cardiovascular system. The focus is just on the simplicity of the overall structure, with a reasonable reflection of the heartfunction. A novel mixed-formulation approach to simulating blood flow in lumped parameters CVS models is outlined that adds minimal complexity, but significantly improves physiological accuracy.

The minimal model is shown to match a Wiggers’ diagram and was also verified to simulate different heartdiseases. The model offers a tool that can be used in conjunction with experimental research to improve understanding of the blood flow.

(6)

(7)

List of Figures

1.1 A schematic picture of the organs in the circulation system . . . 2

1.2 The blood circulation through the heart . . . 3

1.3 An example of a pressure volume diagram together with the ESPVR and the EDPVR lines . . . 5

1.4 The aorta and its hydraulic and electrical representation . . . 7

1.5 a: A modelling lab which consider only the simplest Windkessel method. b: A three elements Windkessel circuit . . . 8

1.6 2-element Windkessel circuit . . . 9

1.7 A Wiggers’ diagram . . . 10

2.1 A simple CVS of a human . . . 14

2.2 A closed loop model of a simple CVS of a human, see figure 2.1 . . . 14

2.3 A possible choice of a cardiac driver function, a = 1, b = 80, c = 0.27, N = 1 . . . 15

2.4 A small part of an artery free of bifurcations . . . 20

2.5 A simple cylindrical artery as a part of the vascular system, where the Γw is the wall of the artery, Γ₁ and Γ₂ are the interfaces with the rest of the system. . . . 23

2.6 Single compartment circuit representation, P pressure, R resistance, C capacitor, Q flow rates . . . 26

2.7 The entire model with in total 12 coupled single compartments . . . 27

2.8 Hydraulic analog of the cardiovascular system. A bifurcation in the systemic circulation is made into a splanchnic and an extrasplanchnic circulation. . . 31

3.1 A simple model of the CVS . . . 41

3.2 A tube free of bifurcations . . . 42

3.3 A cross section of a small artery free of bifurcations . . . 44

3.4 A small part of the cross section . . . 45

3.5 Force balance in the axial direction . . . 46

3.6 An electrical circuit, including a resistor, inductor and capacitor . . . 53

3.7 A simple cardiac driver function, with parameter values: A = 1, B = 80s⁻¹, C = 0.27s and N = 1 . . . 55

3.8 sin² cardiac driver function . . . 56

3.9 sin cardiac driver function . . . 56

4.1 The convergence by different time steps with the Jacobi like method . . . 64

4.2 The convergence for different time steps with the Gauss Seidel like method. In this figure only the 25^th heartbeat is depicted. . . 65

(12)

4.3 The difference in result by using different driverfunctions . . . 68

4.4 Starting with different initial conditions has no influence on the final results . . . 69

4.5 The difference in the results by using inertia . . . 70

5.1 Simulation results from the closed loop model without inertia with our own program 74 5.2 Simulation results from the closed loop model with inertia and ventricular interaction, Results from [SMITH] . . . 75

5.3 A Wiggers’ diagram from a 6 compartment model . . . 76

5.4 A Wiggers’ diagram from a 3 compartment model . . . 78

5.5 Simulating a dystolic disfunction . . . 79

5.6 Effect of diastolic dysfunction causing an increase in ventricle elastance on a PV diagram of the left ventricle [BRWD] . . . 80

5.7 Simulating a systolic disfunction . . . 81

5.8 Effect of systolic dysfunction causing a drop in cintractility on a PV diagram of the left ventricle [BRWD] . . . 81

5.9 Simulating aortic stenosis . . . 82

5.10 A theoretical figure. On the left a normal left ventricle pressure, in the middle a left ventricle pressure caused by aortic stenosis and on the right a left ventricle pressure diagram caused by valvular insufficiency . . . 82

5.11 Simulating valvular insufficiency . . . 83

5.12 Simulating a heart block . . . 84

(13)

Nomenclature

a 1. constant to define the exponential cardiac driver function 2. acceleration

A cross sectional area

b constant to define the exponential cardiac driver function BH compensation term

c constant to define the exponential cardiac driver function C Capacitance

d constant related to the physical properties of the vascular tissues dx longitudinal displacement

dφ infinitesimal angle e(t) cardiac driver function er unit vector in radial direction E Elastance

EW stress

f axisymmetric function F force

h wall thickness

HFB compensation term HFC described threshold I current

l artery length k1 constant

Kr friction parameter

(14)

KHC constant KHG constant L inductance m mass ˆ

p mean pressure over the whole compartment P Pressure

Q 1. intantaneous charge on the capacitor 2. blood flow

Q mean flow rate over the whole districtˆ r internal radius

r₀ reference radius ˆ

r radial direction R Resistance s velocity profile t time

T time interval THF constant u velocity

¯

u mean velocity x axial direction y variable V Voltage

Volume

wp wave propagation β₀ constant

β set of coefficients related to the mechanical and physical properties γ constant

Γw wall of the artery

Γ_1,2 interface with the rest of the system

(15)

ǫ rate of change

η vessel wall displacement θ circumferential coordinate λ constant

µ viscosity

ν 1. kinematic viscosity

2. Poisson ratio of the artery wall ρ density

σ surface stress Φ C¹ function

ψ momentum flux correction coefficient ω heart rate

ωr interaction between fluid and wall

A general axial section P portion of the tube S general axial section V the whole district

CO Cardiac Output

CVS CardioVascular System

EDPVR End Diastolic Pressure Volume relationship ESPVR End Systolic Pressure Volume relationship FE Finite Elements

HR Heart Rate

PRU Peripheral Resistance Unit PV Pressure-Volume

SV Stroke Volume

ZPFV Zero Pressure Filling Volume

(16)

(17)

Chapter 1 Introduction

Death cause number one in the Western world is cardiovascular disease [WE]. Therefor, there is a growing interest in the mathematical and numerical modelling of the human CVS (cardiovascular system). Cited to this interest, much research is devoted to complex three dimensional simulations able to provide sufficient details of the flow field to extract local data such as wall shear stresses. However, these computations are still quite expensive in terms of human re- sources needed to extract the geometry and prepare the computational model and computing time. Since bioengineers and medical researchers do not need the flow in such detail every- where and less detailed models have demonstrated their ability to provide useful information at a reasonable computational cost, further research is done in the description of the CVS in less detailed models. In the less detailed models there must be the possibility to include a small piece of the CVS as a three dimensional model.

In this thesis we do research to a model which 1. Is simple,

2. Needs little computational time and

3. Can accurately reflect a small part of the human CVS.

To create such a model, we start with a literature study to other CVS models. Next the knowledge of others will be used to create a minimal model. Finally, this minimal model will be tested by using the outcomes of models of others and a standard Wiggers’ diagram.

Before starting with the research to different models, we give an introduction in the physiology and in the modelling techniques of blood circulation systems.

1.1 Physiology of the blood circulation system

In the blood circulation system the blood flows through the bloodvessels and is pumped around by the heart. In this section in short there will be an introduction in the physiology of the blood circulation system.

(18)

Figure 1.1: A schematic picture of the organs in the circulation system 1.1.1 The blood circulation system

The blood circulation system consists of a pulmonary- or lungcirculation and a systemic- or bodycirculation. The lungcirculation starts in the heart and provides the lungs of oxygenpoor blood and returns back, with oxygenrich blood, to the heart. Next, the heart pumps the blood into the bodycirculation and provides all the organs of blood before the blood returns to the heart as oxygenpoor blood. One circulation takes about 0.8sec.

1.1.2 The bloodvessels

In the circulation system blood flows through a system of bloodvessels, the vascular system. The bloodvessels are divided into three different groups, the arteries, the capillaries and the veins.

Arteries

After the heart pumps blood away, all of the blood pumped out of the heart (SV - stroke volume) flows into the main artery, the aorta. Most of the SV flows at once into the arterial system to all the organs (see fig 1.1.1). A small part of the SV will be stored in the aorta. Hereby the elastic wall of the aorta will be stretched. When the heart is at rest the aorta contracts and pumps the rest of the SV away. The heart pumps the blood into the arterial system with a pressure of about 120 mmHg (systole). The pressure generated by the aorta is about 80mmHg (diastole).

The blood flows with a velocity of about 4m/sec out of the heart. The velocity in the legs is about 10m/sec. This difference can be explained by the decreasing elasticity of the arteries. Further, the pressure can be influenced by the lumen through narrowing and enlarging (resistance regulation).

(19)

Capillaries

When the blood flows through the organs, the arteries split up in a large system of small arteries, called the capillaries. Because of the large system of small arteries, there is a low presure and a small velocity in the capillaries. The velocity in the cappilaries is about 0.3mm/sec. The advantage of this low pressure and small velocity is that the walls can be thin and the transfer of nutrients with the organs is easy. After the split the capillaries come together in the veins.

Veins

In the veins the flow resistance and pressure drop is small. There are three kinds of mechanisms to pump the blood back to the heart.

• A musclepump For the veins which receive blood from the muscles, there is the muscle pump. By the contraction of the muscle, the vein will be suppressed. By means of a valve the blood in the vein will be pressed in the right direction.

• Arterial-Vein coupling When an arterie and a vein are close to each other the same happens as with a musclepump, but now with an artery which applies pressure on the vein.

• Breath By the underpressure in the chest, the hollow veins are working as a suction-pipe.

1.1.3 Heart

Figure 1.2: The blood circulation through the heart

The heart is a muscle which contains four chambers. By the periodic contraction and relaxation of the muscle, the heart can function as a pump. The four chambers can be split up in two

(20)

atrium-ventricle valve aorta valve

contraction time closed closed

systolic phase

ejection time closed open

cardiaccycle

relaxation time closed closed

diastolic phase

filling time open closed

Table 1.1: Summary of the working of the valves

separate atrium-ventricle couples, a left and the right part. The left atrium-ventricle pair pumps blood through the bodycirculation and the right atrium-ventricle pair pumps blood through the lungcirculation. The left and the right part are working synchron.

Working of the heart

The process of periodic contraction and relaxation of the muscle can be divided in four time periods:

• Contraction time The contraction of the heartmuscle causes a strong increase of the pressure in the ventricle. At this moment, the atrium-ventricle valve and aorta valve are closed.

The volume will be the same (iso-volumetric contraction).

• Ejection time The bloodpressure will be the same as in the artery. The pressure in the ventricle is still increasing. The aorta valve is open. The muscle of the ventricle is con- tracting, the volume will decrease. The valve closes as soon as the blood flows in the wrong direction.

• Relaxation time Relaxation of the heartmuscle. The atrium-ventricle valve is still closed.

There is no change in volume (iso-volumetric relaxation).

• Filling time When the bloodpressure of the ventricle drops beneath the bloodpressure of the atrium, then the atrium-ventricle valve opens.

The contraction time and the ejection time together are called the systolic phase. The relaxation time and the filling time are called the diastolic phase.

1.1.4 Cardiac function

The performance of the heart is indicated by the cardiac function. Three of the most common indicators are the pressure-volume diagram, the cardiac output and preload and afterload. These indicators are used by health professionals to study the patient condition.

The pressure-volume diagram

The pressure-volume (PV) diagram is used to explain the pumping mechanics of the ventricle.

(21)

The two main characteristics of the PV diagram are the lines plotting the End Systolic Pressure- Volume relationship (ESPVR) and the End Diastolic Pressure-Volume relationship (EDPVR), which define the upper and lower limits of the cardiac cycle, respectively.

The cardiac cycle is divided in four parts. In the cardiac cycle, the four time periods of the

Figure 1.3: An example of a pressure volume diagram together with the ESPVR and the EDPVR lines

pumping process can be recognized.

The EDPVR is a measure of the capacitance (C) of the ventricle. The capacitance, defined as the inverse of the elastance (E), is the common term used to describe the PV relationship of an elastic chamber. The ESPVR gives a measure of cardiac contractility, or the strength of contraction, which is defined as the rate at which the heartmuscle reaches peak wall stress.

When there is a diastolic failure, the compliance of the heart wall decreases. Even so, when there is a reduction in contractility the slope of the ESPVR line decreases.

Cardiac Output

The main measure of blood flow on a beat by beat basis is the stroke volume (SV). The SV is defined as the amount of blood pumped from the ventricle during one heart beat. For a more general measure, the cardiac output (CO), is defined as the amount of blood pumped into the aorta, from the left ventricle, in litres per minute. Therefore, the CO is equal to the product of the SV and heart rate (HR):

CO = SV × HR (1.1)

The CO is used to define the capability of the heart to pump nutrient rich blood to the peripheral tissues. The equation for the CO highlights the important dependence on the SV and the HR. While the HR is driven by the sympathetic nervous system, the stroke volume is dependent on the function of the heart muscle as well as on the ventricle preloads and afterloads.

(22)

Preload and Afterload

Preload and afterload are generally intended to be measures of ventricular boundary conditions, indicating the state of the ventricle before and after contraction, respectively. Preload is a measure of the muscle fibre length, immediately prior to contraction, while afterload is a measure of the cardiac muscle stress required to eject blood from a ventricle.

1.2 Cardiovascular System Modelling

Most of the modelling systems for the human CVS can be divided into Finite Element (FE) or Pressure-Volume (PV) approaches. The FE approach involves breaking down parts of the CVS in great detail and utilizing FE calculations to simulate these parts. The PV approach is a simpler method, by grouping parameters and making assumptions to simplify the model as much as possible, while still attempting to simulate the essential dynamics.

1.2.1 Finite Elements Approach

With FE techniques it is possible to get micro-scaled results that can theoretically be very accurate both in trend and magnitude. This kind of approach needs a micro-scale measurement of the mechanical properties such as the elastic properties and dimensions. With this measurements, FE equations can simulate the dynamics of the component being modelled on a micro-scale.

The FE approach on micro-scale is helping to improve understanding. Examples of models using FE techniques are the micro-scale structures of the heart in [NIGRSMHU], [LEHUSM]

and [STHU], or the attempts to model the complex fluid flow dynamics in the heart, particularly around the heart valves in [PEQU] and [GLHUMC]. Although these micro-scale results exists, FE techniques have a lack of flexibility which make them not suitable for patient-specific, rapid diagnostic feedback. Further, it is not feasible to obtain the detailed specific measurements from a living patient, so a model of a specific patient is difficult. Finally, these micro-scale calculations require significant computation time, making these models unsuitable for immediate feedback.

1.2.2 The Pressure-Volume Approach

PV methods are lumped parameter modelling methods where the CVS is divided into a series of elements simulating elastic chambers and blood flow, separately. The elastic chamber elements model the PV relationship in a section of the CVS, such as a ventricle, an atrium, or a peripheral section of the circulation system such as the arteries or veins. All these separate elastic chamber elements are connected by the fluid flow elements which represent blood flow through different parts of the circulation system.

For the modelling of the CVS there exist hydraulic and electrical analogs, see figure 1.4. In the next sections we tell something more about the connection between the electric and hydraulic analog. For the explanation of the hydraulic analogs Windkessel circuits are used. The usage of Windkessel circuits is because most of the PV approaches utilize these circuits.

1.2.3 Windkessel circuit

For most of the calculations the hydraulic formulas are used. To show the connection between the

(23)

Figure 1.4: The aorta and its hydraulic and electrical representation circulation hydraulic electrical

element analog analog

blood flow flow rate current blood pressure pressure voltage pumping function compliance capacitance

vessels viscosity resistance large arteries inertia inductance

Table 1.2: The hydraulic and electric analogs

are circuits which describe the load faced by the heart in pumping blood through the systemic arterial system and the relation between blood pressure and blood flow in the aorta.

One of the first descriptions of a Windkessel circuit was given by the German physiologist Otto Frank in the article ”Die Grundform des Arteriellen Pulses”, published in 1899. In this article Frank compared the heart and systemic arterial system with a closed hydraulic circuit comprised of a waterpump connected to a chamber. The hydraulic circuit is completely filled with water, except for a pocket of air in the chamber. When water is pumped into the system, the water compresses the air in the pocket and pushes water back in the pump. The compressibility of the air in the pocket simulates the elasticity and extensibility of the major artery. This is known as the arterial compliance. The resistance which the water encounters by flowing through the Windkessel and returning back to the pump, simulates the resistance which the blood flow encounters by the blood flowing through the arteries. This process is known as the peripheral resistance.

A Windkessel circuit can consist of a varying number of elements. The simplest model consists of two elements (see figure 1.2.3), namely a compliance and a peripheral resistance. By using the basic laws of an electrical circuit (Ohm’s law and Kirchhof’s laws), the Windkessel

(24)

Figure 1.5: a: A modelling lab which consider only the simplest Windkessel method. b: A three elements Windkessel circuit

model can be described by a mathematical model.

According to Ohm’s law, the drop in electrical potential across the resistor is IRR and the drop in electrical potential across the capacitor is Q/C, where Q is the instantaneous charge on the capacitor and ^dQ_dt = IC. From Kirchhof’s voltage law, the net change in electrical potential around each loop of the circuit is zero; therefor V (t) = IRR and V (t) = Q/C. From Kirchhof’s current law, the sum of currents into a junction must equal the sum of currents out of the same junction: I(t) = IC+ IR. By now, the current in the capacitor is given by IC = C(dV /dt). If we now substitute IC and IRfrom above into Kirchhof’s current law then we finally get an electric mathematical model which describes the 2-element Windkessel model:

I(t) = CdV (t)

dt +V (t)

R (1.2)

In terms of the physiological system, I(t) is the blood flow from the heart to the aorta, V (t) is the blood pressure in the aorta, C is the arterial compliance and R is the peripheral resistance in the arterial system. In physiological terms the hydraulic mathematical model reads:

Q(t) = CdP (t)

dt +P (t)

R (1.3)

Now, we use the hydraulic equivalent to evaluate what happens during diastole. During diastole, there is no inflow, so Q(t) = 0 and an exact solution exists:

P (t) = P (0)e^−RCt (1.4)

(25)

Figure 1.6: 2-element Windkessel circuit

A modification of the 2-element Windkessel circuit is obtained by including an inductor in the main branch of the circuit, as can be seen in figure (1.5b). This inductor simulates inertia of the fluid in the hydrodynamic model. The mathematical model of this 3-element Windkessel circuit can be found by using that the drop in electrical potential across an inductor with inductance equals VL= L(dIL(t)/dt) and Kirchhof’s law, IL= IR+ IC:

IL(t) = IR(t) + IC(t) = ^V_R^R + C^dV_dt^C = ^V^(t)−V_R ^L + C^d^{(V (t)−V}_dt ^L⁾ = ^V_R^(t) −^VR^L + C^dV_dt^(t) − C^dVdt^L

= ^V_R^(t) −^LR dI(t)

dt + C^dV_dt^(t) − CL^d²_dt^I^L²^(t)

⇔

I_L(t) +_R^L^dI^L_dt^(t) + CL^d²_dt^I^L₂^(t) = ^V_R^(t) + C^dV_dt^(t) In this case, the hydraulic equivalent reads:

Q(t) + L R

dQ

dt + LCp

d²Q

dt² = P (t) R + Cp

dP (t)

dt (1.5)

1.2.4 Wiggers’ diagram

The Wiggers’ diagram depict the pressure and volume in the heart and the ejecting activity of the heart:

• The first diagram shows the electrocardiogram (ecg). We will not use this ecg.

(26)

Figure 1.7: A Wiggers’ diagram

• The second diagram represents the pressure in the left atrium, left ventricle and in the aorta.

The top at a of the left atrium and ventricle is the result of the contraction of the atrium.

Next the ventricle contraction occurs. During ventricle ejection, the atrium is pulled. The result of this pulling is a pressure drop in the atrium, showed at the x-top. Right before the x-top, a c-top in the pressure of the left atrium is denoted. This c-top is caused by the opening of the aortic valve. In the second part of the ejection of the left ventricle, the pressure in the atrium is increasing due to the filling with blood until the mitral valve is opened for a fast filling of the ventricle during the isometric relaxation. This opening causes the y-top.

• The third diagram depicts the volume in the left ventricle and the flow velocity in the aorta.

(27)

• The fourth diagram is similar to the second diagram, with the difference that this diagram depicts the values for the right atrium and ventricle. As can be seen, the pressure in the right side of the heart is lower then in the left side of the heart.

• In the last diagram a phonocardiogram is showed. The first (I) noise reflects the closing of the mitral valve, the second (II) noise reflects the closure of the aortic valve. The first and second noise give exactly the duration of the relaxation and contraction. The systolic phase starts at the beginning of the first noise till the beginning of the second noise. The diastolic phase starts at the beginning of the second noise till the beginning of the first noice.

What is the connection with our reseach? As said before, we want to use the Wiggers’ diagram as reference material for a specific person. The phonocardiogram can be used to measure the time needed for the different heart periods. The other three graphs can be used as reference for our own model.

1.3 Summary

In this section the main goals of this report have been outlined, namely creating a human CVS model on a low level. Further it must be possible to describe a small part of the CVS detailed.

To make it possible to create a good model for the human CVS, we gave an introduction in the physiology of the circulation system. Finally, we gave an introduction about the possibilities of CVS modelling. We introduced two approaches, a FE approach and a PV approach, together with their advantages and disadvantages.

(28)

(29)

Chapter 2 Literature study

In this literature study we are going to search for existing models. The models found in literature, which at a first glance satisfy the requirements mentioned in the introduction, will be described.

Based on the summaries we will make decisions about the model we will use in further studies.

Every section contains the summary of one model. The conclusions will be presented in the final section.

2.1 Minimal Heamodynamic Modelling of the Heart & Circu- lation for Clinical Application [SMITH]

2.1.1 General description

This is the title of the thesis that is presented by Bram W. Smith for the degree of Doctor of Philosophy in Mechanical Engineering at the University of Canterburry, Christchurch, New Zealand. In this thesis Smith has the intention to make a model which contains:

• A closed-loop, stable model with minimal complexity and physiologically realistic inertia and valve effects.

• A model parameters that can be relatively easily determined or approximated for a specific patient using standard, commonly used techniques.

• A model that can be run on a standard desktop computer in reasonable time (eg. in the order of 1-5 minutes)

• Accurate prediction of trends

The model presented a hydraulic, 0D, 6 compartment model, which intends to simulate the essential haemodynamics of the CVS including the heart, and the pulmonary and systemic circulation systems. Figure 2.1 shows a simplified diagram of the human circulation system with in the middle the human heart. Figure 2.2 presents a closed-loop model of the same human CVS. As can be seen in figure 2.2, the closed-loop model contains compartments which are connected by resistors and inductors in series and can be seen as a Windkessel circuit.

For the pulsation of the heart, Smith uses a cardiac driver function e(t). This cardiac driver function utilizes the ESPVR and EDPVR (see 1.1.4) as the upper and lower limits of cardiac

(30)

Figure 2.1: A simple CVS of a human

Figure 2.2: A closed loop model of a simple CVS of a human, see figure 2.1

chamber elastance. The profile of the driver function represents the variance of elastance between minimum and maximum values during a single heart beat:

e(t) =

N

X

i=1

a_ie^−bⁱ^(t−cⁱ⁾², (2.1)

where the ai, bi, ci and N are parameters that determine the shape of the driver profile. For his simple model he takes a = 1, b = 80s⁻¹, c = 0.27s and N = 1. See for the shape figure 2.3.

With respect to the heart Smith makes some assumptions, which will be described below.

2.1.2 Assumptions

The first assumption that Smith makes, is that blood, which flows through the CVS, is approximated as flow through a tube. The flow rate equations have directly been derived from the

(31)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Time t (sec)

e(t)

profile of the cardiac driver function

Figure 2.3: A possible choice of a cardiac driver function, a = 1, b = 80, c = 0.27, N = 1 Navier-Stokes equations in cylindrical coordinates:

∂ux

∂t + ur

∂ux

∂r +uθ

r

∂ux

∂θ + ux

∂ux

∂x = −1 ρ

∂P

∂x + ν 1 r

∂

∂r

r∂ux

∂r

+ 1

r²

∂²ux

∂θ² +∂²ux

∂x²

(2.2) where ux, ur and uθ are the longitudinal, radial and angular velocities, respectively, P is the modified pressure relative to hydrostatic, ρ is the density and ν is the kinematic viscosity.

The next assumptions are standard assumptions, which will be applied to all equations governing fluid flow.

• Blood is assumed to be incompressible, so ρ is constant.

• For the heart, the fluid is assumed to behave in a continuous, Newtonian manner with constant viscosity (µ is constant).

• The arteries are assumed to be rigid with a constant cross sectional area (^∂r∂x = 0). This assumption fits with standard Windkessel circuit design involving a rigid pipe and an elastic compartment in series. The rigid tube simulates the fluid dynamics, while the elastic compartment simulates the compliance of the artery walls.

• Laminar uni-directional axi-symmetric flow is assumed (u^r = 0, uθ = 0 and ^∂u_∂θ^x = 0).

Although turbulence can occur around the valves, it takes time to develop, and is assumed not to affect the flow profile significantly.

• The flow is assumed to be fully developed along the length of the tube meaning the velocity profile is constant with respect to x (^∂u_∂x^x = 0).

• Pressure is assumed constant across the cross-sectional area and the pressure gradient is constant along the length of each section so that the pressure gradient is a function of time only (^∂P_∂x(t)).

With these assumptions equation (2.2) reduces to the following equation ρ∂u(r, t)

∂t = −∂P

∂x(t) +µ r

∂

∂r

r∂u(r, t)

∂r

, (2.3)

(32)

where µ is the viscosity (µ = νρ) and u(r, t) is the velocity in the x-direction (ux(r, t)) as a function of radius and time only.

2.1.3 Mathematical Model

For the description of the compartments Smith has the choice between two different equations.

Both equations are derived from equation (2.3). The first equation used is Poiseuille’s equation for flow rate, assuming constant resistance and no inertial effects:

Q(t) = P₁(t) − P2(t)

R , (2.4)

where R = _πr^8µl4

0 is the resistance. The second equation includes inertial effects and constant resistance

LdQ

dt = P₁− P2− QR, (2.5)

where L = _πr^ρl2

0 is inertia and R = _πr^8µl4

0 constant resistance. Around the heart where are big flow differences the second equation will be used. Far from the heart flows a nearly constant flow, so the first equation will be used.

2.1.4 Parameters

Smith uses the following parameters for his tests.

(33)

Description Symbol Value Blood properties

Blood Density ρ 1050kg/m³

Blood Viscosity µ 0.004N S/m²

Blood Kinematic Viscosity ν 3.8 10⁻⁶m²/s Stressed Volume of blood in CVS Vtot 1500ml Unstressed Volume of blood in CVS 4000ml

Artery properties

Internal Artery Radius r₀ 0.0125m

Artery length l 0.2m

Compartment properties

Chamber Elastance E_es 1N/m⁵

EDPVR Volume Cross-over V₀ 0m³

ESPVR Volume Cross-over Vd 0m³

Constant λ 23000m⁻³

Heart Rate ω 1.33beats/sec

Constant a 15N/m²

Table 2.1: Constants used in a single compartment simulation

PARAMETERS E_es V_d V₀ λ P₀

Units 10⁶N/m⁵ 10⁻⁶m³ 10⁻⁶m³ m⁻³ N/m²

Left Ventricle free wall (lvf) 100 0 0 33000 10

Right ventricle free wall (rvf) 54 0 0 23000 10

Septum free wall (spt) 6500 2 2 435000 148

Pericardium (pcd) - - 200 30000 66.7

Vena-cava (vc) 1.3 0 - - -

Pulmonary Artery (pa) 72 0 - - -

Pulmonary Vein (pu) 1.9 0 - - -

Aorta (ao) 98 0 - - -

Table 2.2: Mechanical properties of the heart and circulation system

(34)

Resistance Inertance

Parameter N s/m⁵ N s²/m⁵

Mitral Valve (mt) 6.1 10⁶ 1.3 10⁴

Aortic Valve (av) 2.75 10⁶ 5 10⁴

Tricuspid Valve (tc) 1 10⁶ 1.3 10⁴

Pulmonary Valve (pv) 1 10⁶ 2 10⁴

Pulmonary Circulation System (pul) 9.4 10⁶ N/A Systemic Circulation System (sys) 170 10⁶ N/A

Table 2.3: Hydraulic properties for flow between compartments

Description Symbol

Elastance of Vena-cava E_vc = 1.29 10⁶N/m⁵ Elastance of Left Ventricle P_0,lvf = 9.07N/m² Elastance of Pulmonary Artery Epa = 44.5 10⁶N/m⁵ Elastance of Pulmonary Vein E_pu= 0.85 10⁶N/m⁵ Elastance of Right Ventricle P_0,rvf = 20.7N/m² Elastance of Aorta Eao = 98 10⁶N/m⁵ Resistance of Tricuspid Valve Rtc= 3.3 10⁶N s²/m⁵ Resistance of Pulmonary Valve R_pv= 1 10⁶N s²/m⁵ Resistance of Pulmonary Circulation Rpul= 19.3 10⁶N s²/m⁵ Resistance of Mitral Valve Rmt= 2.33 10⁶N s²/m⁵ Resistance of Aortic Valve Rav= 5.33 10⁶N s²/m⁵ Resistance of Systemic Circulation R_sys = 139.6 10⁶N s²/m⁵ Contractility of Left Ventricle Ees,lvf = 377 10⁶N/m⁵ Contractility of Right Ventricle Ees,rvf = 87.8 10⁶N/m⁵

Table 2.4: Parameter values for the closed loop model

(35)

2.1.5 Conclusions

Smith performs different tests to verify all the specific possibilities of his model. One of the tests is the comparison with a Wiggers’ diagram (for explanation of a Wiggers’ diagram see paragraph 1.2.4). In table 2.5, Smith compares his results with the Wiggers’ diagram. After

Model Target

Variable Value Value % Error Pressure in Aorta

Amp Pao 40mmHg 41.407 3.5%

Avg Pao 100mmHg 119.168 19.2%

Pressure in Pulmonary Artery

Amp Ppa 17mmHg 20.414 20.1%

Avg Ppa 16.5mmHg 20.314 23.1%

Volume in Left Ventricle

Amp Vlv 70ml 69.508 −0.7%

Avg Vlv 80ml 84.042 5.1%

Volume in Right Ventricle

Amp Vrv 70ml 69.569 −0.6%

Avg Rrv 80ml 121.185 51.5%

Pressure in Pulmonary Vein

Avg Ppu 2mmHg 10.112 405.6%

Pressure in Vena-cava

Avg Pvc 2mmHg 1.050 −47.5%

Table 2.5: Comparison of the results with the Wiggers’ diagram

doing different tests, Smith comes with the following conclusions:

• The blood flow rate is primarily dependent on the pressure gradient across the resistor.

If the effects of inertia are either ignored or negligible, the equation for flow rate can be calculated using Poiseuilles equation (2.4). Poiseuilles equation assumes incompressible, Newtonian, laminar, axi-symmetric, fully developed flow through a rigid tube of constant cross-section.

• Tests prove the stability of the closed loop CVS model.

• The model is seen to capture the major dynamics of the CVS including the variations in left ventricle pressure, aortic pressure and ventricle volume.

• The decrease in cardiac output is in good agreement with readily available clinical data.

• The results show the capability of the presented approach to create patient specific models.

(36)

2.2 Reduced and multiscale models for the human cardiovascu- lar system; one dimensional model [FORVEN]

Formaggia and Veneziani wrote this report as collection of the notes of the two lectures given by Formaggia at the 7^th VKI Lecture Series on ”Biological fluid dynamics” held on the Von Karman Institute, Belgium, on May 2003. They give a summary of some aspects of the research aimed at providing mathematical models and numerical techniques for the simulation of the human CVS.

At first they derive an one dimensional model. Hereto, they start with the mathematical

Figure 2.4: A small part of an artery free of bifurcations

description of a small part of an artery free of bifurcations. They assume that the small part of the artery can be described by a straight cylinder with a circular cross section. For the description of the flow through this straight cylinder the Navier-Stokes equations are used and integrated over a generic cross section. Starting parameters are the time interval T = (0, t₁) and the vessel length x = (0, l).

2.2.2 Assumptions

Describing a small part of the artery as a straight cylinder with the Navier-stokes equations is too expensive, so some simplifying assumptions are made:

1. All quantities are independent of the angular coordinate θ. As a consequence, every axial section x = const remains circular during the wall motion. The tube radius r is a function of x and t.

2. The wall displaces along the radial direction solely, thus at each point at the tube surface they may write η= ηer, where η = r − r0 is the displacement with respect to the reference radius r0.

3. The vessel will expand and contract around its axis, which is fixed in time. This hypothe- sis is indeed consistent with that of axial symmetry. However, it precludes the possibility

(37)

of accounting for the effects of displacements of the artery axis such as occuring in the coronaries because of the heart movement.

4. The pressure P is constant on each section, so it only depends on x and t.

5. The body forces are neglected.

6. The velocity components orthogonal to the x axis are negligible compared to the component along x. The latter is indicated by ux and its expression in cylindrical coordinates is supposed to be of the form

ux(t, ˆr, x) = ¯u(t, x)s(ˆrr⁻¹(x)) (2.6) where ¯u is the mean velocity on each axial section and s : R → R is a velocity profile, which must be chosen such that R₁

0 s(y)ydy = ¹₂. The fact that the velocity profile does not vary in time and space is in contrast with experimental observations and numerical results carried out with full scale models. However, it is a necessary assumption for the derivation of the reduced model. One may then think of s as being a profile representative of an average flow configuration.

Finally, a momentum-flux correction coefficient is defined by:

ψ = R

Su²_xdσ A¯u² =

R

Ss²dσ

A .

where A is the cross sectional area and S the general axial section.

2.2.3 Mathematical model

With all these assumptions, the main variables are:

• Q the mean flow, defined as

Q = Z

S

uxdσ = A¯u;

• A the surface area of an axial section;

• P the pressure.

When ψ is taken constant, the reduced model looks like ( _∂A

∂t +^∂Q_∂x = 0

∂Q

∂t + ψ_∂x^∂ (^Q_A²) +^A_ρ^∂P_∂x + K_r(^Q_A) = 0 (2.7) for x ∈ (0, l), t ∈ T , where K^r = −2πνs^′ is a friction parameter and s^′ the derivative of the velocity profile.

(38)

parameter

input pressure amplitude 20 × 10³dyne/cm²

FLUID viscosity, ν 0.035poise

Density, ρ 1.021kg/m³

STRUCTURE Wall Thickness, h 0.05cm Reference Radius, r₀ 0.5cm

Table 2.6: Parameters used in the one dimensional model 2.2.4 Parameters

For the model specific parameters see [FORVEN, page 1.37].

2.2.5 Conclusions

Before a test can be done, a velocity profile has to be chosen. Hereto several options exist.

Formaggia and Veneziani chose for the parabolic profile s(y) = 2(1−y²). This profile corresponds with the Poisseuille solution characteristic of steady flow in circular tubes. The parabolic profile is a variant of the profile most used: s(y) = γ⁻¹(γ + 2)(1 − y^γ). This power law profile is most used, since it has been found experimentally that the velocity profile is, on average, rather flat.

Now that a velocity profile has been selected, three sets of tests are distinguished. The first series of tests are focussed on the single artery, the second series of tests are done with a coupling of 55 main arteries and the last series of tests are an improvement of the second series of tests.

The improvements in the third test are made by taking inertia of the wall into account. For the results see [FORVEN].

Comparing the results with literature, Formaggia and Veneziani conclude that there is little agreement with reality. This can be explained by the chosen model. The model is namely formed by a closed network with a high-level of inter-dependency. In this model the flow dynamics of the blood in a specific vascular district is stricly related to the global, systemic dynamics.

However, in [ARFEL], it is shown that even a strong reduction in the vascular lumen in a carotid bifurcation does not mean a relevant reduction of the blood supply to the brain.

So, to make the results more realistic another way of coupling parts of a high inter-dependence model must be found. The next section describes how Formaggia and Veneziani make use of a Windkesselcircuit to accomplish this.

2.3 Reduced and multiscale models for the human cardiovas- cular system;lumped parameters for a cylindrical compliant vessel [FORVEN]

In the previous section Formaggia and Veneziani developed a 1D model. After doing different tests they conclude that the model has a high level of interdependency and that there is little agreement reflecting reality. So, they continu their research aimed at providing mathematical models and numerical techniques for the simulation of the human CVS, by focussing on coupling techniques. In this section Formaggia and Veneziani are describing a mathematical model of the CVS which couples a local system with a systemic model. The local system is based on the

(39)

solutions of the incompressible Navier-Stokes equations possibly coupled with the dynamics of the vessel wall, while the systemic model is based on a one-dimensional system or on a lumped parameters model. The lumped parameters model is based on the solution of a system of ordi- nary differential equations for the average mass flow and pressure.

For the systemic model, a choice can be made between a one dimensional model and a lumped parameters model. As one dimensional model you can think of a model like the one described in section 2.2. A lumped parameters model is described below.

A lumped parameter models

A lumped parameters model provides a systemic description of the main phenomena related to the circulation at a low computational cost. An effective description of this model is by dealing with separate ’compartments’ and their interaction. To develop a lumped parameters model

Figure 2.5: A simple cylindrical artery as a part of the vascular system, where the Γ_w is the wall of the artery, Γ₁ and Γ₂ are the interfaces with the rest of the system.

they start with a seperate compartment, for which they consider a simple cylindrical artery, see figure 2.5. In this circular cylindrical domain, the axial section equals A(t, x) = πr²(t, x) where r(t, x) is the radius of the section at x. With this consideration they want to form a simplified model, therefor some assumptions have been introduced, as described next.

2.3.2 Assumptions

As in the one dimensional model of section 2.2, Formaggia and Veneziani start with the Navier- Stokes equations and make the same standard assumptions. After these standard assumptions a 1D model is left.

( _∂A

∂t +^∂Q_∂x = 0

∂Q

∂t + ψ_∂x^∂

Q² A

+^A_ρ^∂P_∂x + KrQ

A = 0 (2.8)

In order to close the system, a further equation, which is provided by the constitutive law for the vessel tissues, is needed. So another assumption is made:

• The vessel wall displacement η is related to the pressure P by an algebraic linear law. By following [FOVE], they take:

(P − P^ext) = d(r − r⁰) = β0

√A −√ A₀

A₀ , (2.9)

where Pext is a constant reference pressure, A₀ = πr²₀ a constant reference area, d is a constant related to the physical properties of the vascular tissues and β₀ = A₀d/√

π.

(40)

Now, the one dimensional model (2.8) will be integrated along x ∈ (0, l):







k₁l^dˆ_dt^p + Q₂− Q1= 0 l^{d ˆ}_dt^Q+ ψh_Q₂

2

A2 −^Q_A²¹₁i +Rl

0

hA ρ

∂P

∂x + KRQ A

i

dx = 0 (2.10)

In this system, ˆp is the mean pressure over the whole compartment and k₁ a constant. Since this is not a satisfactory model, because it is not linear, some more assumptions are needed.

• The quantity _Q2 2

A2 −^Q_A₁²¹

is so small compared to the other terms in short pipes that it can be discarded.

• The variation of A with respect to x is small compared to that of P and Q, so the integral in 2.10 will be approximated

Z l 0

A ρ

∂P

∂x + K_RQ A

dx ≈

Z l 0

A0

ρ

∂P

∂x + K_RQ A0

dx

2.3.3 Mathematical model

With all these assumptions Formaggia and Veneziani end op with a system for the lumped parameters description of the blood flow in the compliant cylindrical vessel. It involves the mean values of the flow rate and the pressure over the domain, as well as the upstream and downstream flow rate and pressure values:







k₁l^dˆ_dt^p + Q₂− Q1 = 0

ρl A0

d ˆQ

dt +^ρK_A^R2^l 0

Q + Pˆ 2− P¹ = 0 (2.11)

The final system (2.11) will be represented by a hydraulic analog. In this analog three parameters are used, namely:

• R the resistance induced by the blood viscosity is represented by R = ^ρK_A²^R^l

0 . Assuming a parabolic velocity profile gives

R = 8µl

πr⁴₀ (2.12)

• L the inductance of the flow represents the inertial term in the momentum conservation law and is given by

L = ρl

πr₀² (2.13)

• C the capacitance of the vessel represents the coefficient of the mass storage term in the mass conservation law and is given by

C = 3πr₀³l

2Eh (2.14)

With this notation equation (2.11) becomes

( C^dˆ_dt^p + Q₂− Q1 = 0 L^{d ˆ}^Q + R ˆQ + P − P = 0

(2.15)

(41)

2.3.4 Conclusions

In an analytical test case a completely lumped parameters description of the circulation, providing a reference solution to the systemic level, is given. The aim of this test is to model blood flow behaviour in Ω by the Navier-Stokes equations coupled with the lumped description of the remaining network. Formaggia and Veneziani expect from this test that the presence of a local accurate submodel does not have to modify significantly the results at the systemic level. This is exactly what they obtain numerically. The heterogeneous model is able to compute accurately the velocity and pressure fields in the domain of interest.

In a test case where a 3D-1D coupling is made, it could be concluded that the coupling between a 3D fluid structure model and a 1D reduced model is an effective way to greatly reduce numerical reflections of the pressure waves.

In a last test of clinical interest, the methodology was in particular applied to a recon- structive procedure, the systemic-to-pulmonary shunt, used in cardiovascular paediatric surgery to treat a group of complex congential malformations. The 3D-model includes the shunt, the innominate artery (through which blood flows in) and the pulmonary, carotid and subclavian arteries (through which blood flows out). In this case the lumped model is composed by different blocks describing the rest of the pulmonary circulation, the upper and lower body, the aorta, the coronary system and the heart. This application to the systemic-to-pulmonary shunt, gives a clear idea of what can be obtained using the multiscale methodology.

2.4 Computational modeling of cardiovascular response to or- thostatic stress [HSKM]

The objective of this study is to develop a model of the cardiovascular system capable of simulating the short term (< 5min) transient and steady state heamodynamic responses to head-up tilt and lower body negative pressure. A subobjective of this study is to develop and test a general 0D, 12 compartment model of a CVS that contains the essential features associated with the effects of gravity. The development of the model is not completely their own, but they use the knowledge and formulas of other investigators.

2.4.2 Mathematical Model

The model of [WHFICR] and [DAMA] is based on a closed-loop lumped parameters heamodynamic model with local blood flow to major peripheral circulatory branches. This heamodynam- ical model is mathematically formulated in terms of a hydraulic analog model in which inertial effects are neglected. A single compartment circuit representation which has been used is given in figure 2.6. The equations read

Q₁= Pn−1− Pⁿ Rn

(2.16) Q₂= Pn− Pn+1

R_n+1 (2.17)

Q₃= d

dt[Cn× (Pⁿ− Pbias)] (2.18)

(42)

Figure 2.6: Single compartment circuit representation, P pressure, R resistance, C capacitor, Q flow rates

The flow Q₁ at node Pn splits up like Q₁ = Q₂+ Q₃. Combining these expressions for the flow rates leads to

d

dtPn= Pn+1− Pⁿ

CnR_n+1 +Pn−1− Pⁿ CnRn

+Pbias− Pⁿ

Cn · dCn

dt

+ d

dtPbias (2.19) In total 12 such first order differential equations are used to describe the entire model, see figure 2.7. To solve this model, the authors used a fourth order Runge-Kutta integration routine. The pumping action of the heart is realized by varying the right and the left ventricular elastances according to a predefined function of time (Er and El)

e(t) =











E_dias+^E^sys^−E₂ ^dias

1 − cos

π ^t

0.3√

T(n−1)

0 ≤ t ≤ Ts

Edias+^E^sys^−E₂ ^dias

1 + cos

2π^t^−0.3

√T(n−1) 0.3√

T(n−1)

Ts< t ≤ ³₂Ts

E_dias ³₂T_s< t ≤ T (n)

(2.20)

In this equation Edias and Esys represent the end-diastolic and end-systolic elastance values, respectively. Further T (i) denotes the cardiac cycle length of the ith beat and t denotes the time measured with resect to the onset of ventricular contraction. The systolic time interval, Ts, is determined by the Bazett formula, Ts(n) ≈ 0.3pT (n − 1). The atria are not represented, because their function is partially absorbed into the function of adjacent compartments.

For the change of volume in the compartments the authors refer to experimental observations of [LUD]. In accordance to these observations they model the functional form of the pressure- volume relationships of the venous compartments of the legs, the splanchnic circulation and the abdominal venous compartment with

∆V = 2∆V

π arctan

πC₀ 2∆Vmax

∆Ptrans

, (2.21)

where ∆V represents the change in compartment volume due to a change in transmural pressure ∆Ptrans, ∆Vmax is the maximal change in compartment volume and C₀ represents the

(43)

Figure 2.7: The entire model with in total 12 coupled single compartments

is modified as a function of time to simulate fluid sequestration into the interstitium during orthastatic stress.

2.4.3 Parameters

The parameters given in this section are the resistance, volume and capacitance values for the 12 different compartments.

In this tables the writers make use of the P RU = mmHg.s/ml (peripheral resistance unit) and ZP F V (zero pressure filling volume). Further all values compound with a 71 − 75kg normal male subject and a body surface area of 1.7 − 2.1m² with a total blood volume of 5700ml.

(44)

Resistance values [P RU ]

Rlo R_up1 R_kid1 R_sp1 R_ll1 R_up2 R_kid2 R_sp2

0.006 3.9 4.1 3.0 3.6 0.23 0.3 0.18

R_ll2 Rsup Rab Rinf Rro Rp Rpv

0.3 0.06 0.01 0.015 0.003 0.08 0.01

Table 2.7: Resistance values

ZPFV Capacitance

Compartment ml ml/mmHg

Right ventricle 50 1.2-20 Pulmonary arteries 90 4.3

Pulmonary veins 490 8.4 Left ventricle 50 0.4-10 Systemic arteries 715 2.0

Systemic veins

Upper body 650 8

Kidney 150 15

Splanchnic 1300 55

Lower limbs 350 19

Abdominal veins 250 25

Inferior vena cava 75 2 Superior vena cava 10 15

Table 2.8: Volume and capacitance values

2.4.4 Conclusions

After several tests have been done, it shows that all major heamodynamic parameters generated by the model are within the range of what is considered as physiologically normal in the general population. Representative simulated pressure waveforms are made, too. The conclusions the authors draw after these tests are

1. They assume that the dynamics of the system can be simulated by restricting their analysis to relatively few representative points within the CVS. Although this approach is incapable of simulating pulse wave propagation, it does reproduce realistic values of beat-to-beat heamodynamic parameters.

2. One potential limitation of the heamodynamic system in its present form might be the lack of atria, which are thought to contribute significantly to ventricular filling at high heart rates.

3. The model generates steady state and transient heamodynamic responses that compare well to population-averaged and individual subject data.

modelling of the circulation

Minimal haemodynamic

modelling of the circulation

P.C.I. Spelde

Master Thesis in Applied Mathematics

April 2008

Minimal haemodynamic

modelling of the circulation

P.C.I. Spelde

First supervisors: A.E.P. Veldman and G. Rozema Second supervisor: A.J. van der Schaft

External supervisor: N.M. Naurits

Institute of Mathematics and Computing Science P.O. Box 407

9700 AK Groningen

The Netherlands

Abstract

Contents

List of Figures

Nomenclature

Chapter 1

Introduction

1.1 Physiology of the blood circulation system

1.2 Cardiovascular System Modelling

1.3 Summary

Chapter 2

Literature study

2.1 Minimal Heamodynamic Modelling of the Heart & Circu- lation for Clinical Application [SMITH]

2.2 Reduced and multiscale models for the human cardiovascu- lar system; one dimensional model [FORVEN]

2.3 Reduced and multiscale models for the human cardiovas- cular system;lumped parameters for a cylindrical compliant vessel [FORVEN]

2.4 Computational modeling of cardiovascular response to or- thostatic stress [HSKM]