• No results found

Parametric Study of a Mitral Valve Model for Blood Flow Simulation in the Left Ventricle of the Human

N/A
N/A
Protected

Academic year: 2021

Share "Parametric Study of a Mitral Valve Model for Blood Flow Simulation in the Left Ventricle of the Human "

Copied!
30
0
0

Loading.... (view fulltext now)

Full text

(1)

IN

DEGREE PROJECT TECHNOLOGY, FIRST CYCLE, 15 CREDITS

STOCKHOLM SWEDEN 2020,

Parametric Study of a Mitral Valve Model for Blood Flow Simulation in the Left Ventricle of the Human

Heart

HANNA KYLHAMMAR

ADAM SKÄRBO JONSSON

KTH ROYAL INSTITUTE OF TECHNOLOGY

(2)

Abstract

Technical solutions for diagnosis and treatment of heart diseases is of great interest due to it being one of the leading causes of death in the world. This study focuses on the valve between the left atrium and ventricle, the mitral valve, and the blood flow in the left ventricle of the human heart. A paramet- ric model for patient-specific simulations of mitral valve dynamics developed by Hoffman J. et al. at KTH and Lucor D. at Universit´e Paris-Saclay is used and further studied by focusing on the heart disease mitral valve stenosis.

Solutions are produced with FEniCS-HPC HeartSolver using computational fluid dynamics. Our results show that the more severe the disease is, the higher velocities and vorticity in the ventricle, and the lower the pressure.

Further studies of interest include non-Dirichlet boundary conditions by the mitral valve, to allow mitral regurgitation that was observed in this study.

Sammanfattning

Tekniska l¨osningar f¨or diagnostisering och behandling av hj¨artsjukdomar ¨ar av stort intresse f¨or samh¨allet med anledning att det ¨ar en av de fr¨amsta d¨odsor- sakerna v¨arlden ¨over. Denna avhandling fokuserar p˚a klaffen i hj¨artat mel- lan v¨anstra f¨ormaken och kammaren, mitralisklaffen, och blodfl¨odet i v¨anstra kammaren. En parametermodell f¨or patientspecifika simuleringar, framtagen av Hoffman J. et al p˚a KTH och Lucor D. p˚a Universit´e Paris-Saclay, anv¨ands och utvecklas vidare genom att till¨ampa modellen p˚a hj¨artsjukdomen mi- tralisstenos. L¨osningarna ¨ar baserade p˚a ber¨akningsstr¨omningsdynamik och framtagna med FEniCS-HPC HeartSolver. V˚ara resultat visar p˚a att ju al- lvarligare symptom desto h¨ogre hastighet och vorticitet i kammaren, samt desto l¨agre tryck. Vi rekommenderar f¨orb¨attringar av modellen genom att

¨

andra randvillkoren fr˚an Dirichletvillkor f¨or att till˚ata regurgitation genom mitralisklaffen, d˚a v˚ara resultat pekade p˚a att det var ett m¨ojligt utfall.

(3)

CONTENTS CONTENTS

Contents

1 Introduction 5

2 Background 5

2.1 The Human Heart . . . 5

2.2 The Mitral Valve . . . 5

2.2.1 Mitral Valve Diseases . . . 6

2.3 Numerical Hemodynamic Modeling . . . 6

2.3.1 Navier-Stokes Equations for Moving Boundaries . . . 6

2.3.2 Vorticity and Streamlines . . . 7

2.3.3 Finite Element Method . . . 7

2.4 Curve Fitting . . . 8

2.4.1 Cubic Splines . . . 9

2.4.2 Curve Fitting with Splines . . . 9

2.5 Simulation of a Healthy Mitral Valve . . . 10

3 Method 11 3.1 Mitral Valve Design . . . 11

3.1.1 Model for Mitral Stenosis . . . 16

3.2 Variation of Parameters . . . 17

3.3 Patient Specific Data . . . 17

3.4 Software . . . 19

3.4.1 FEniCS-HPC HeartSolver . . . 19

3.4.2 DOLFIN-HPC . . . 20

3.4.3 DOLFIN-post and ParaView . . . 21

4 Results 21 4.1 Velocity . . . 21

4.2 Pressure . . . 22

4.3 Vorticity . . . 23

4.4 Streamlines . . . 25

5 Discussion 26 5.1 Conclusion . . . 26

5.2 Choice of Diseases to Study . . . 27

5.3 Ethical Aspects . . . 27

5.4 Further Studies . . . 27

5.4.1 Shortening Axes with a Constant . . . 27

5.4.2 Negative Flow . . . 28

Appendices 30

(4)

CONTENTS CONTENTS

A Ramping Coefficients 30

(5)

2 BACKGROUND

1 Introduction

Proposed in previous work[1, 2] is a parametric model of the mitral valve (MV).

The aim of this paper is to perform a parametric analysis of an existing MV model and the disease mitral stenosis in a simulated environment. Simulation of the dy- namics of the blood flow in the left ventricle uses images of the ventricle generated from echocardiography that provide needed moving boundaries, volume data and sophisticated model of the mitral valve [3]. This type of fluid dynamics simulation requires high computing power. Therefore, the simulations are made with FEniCS- HPC HeartSolver [1, 4] running on the supercomputing cluster Beskow provided by PDC KTH [5].

2 Background

In this section we will provide the reader with a basic understanding of the heart and mitral valve, as well as an introduction to the equations the simulations are based on and how they are solved in this particular problem.

2.1 The Human Heart

Blood is oxygenated by the lungs and pumped through the body by the heart. The human heart consists of two atria and two ventricles, left and right. The path of the blood through the heart is: from the lungs to the left atrium to the left ventricle to the aorta, then from the body to the right atrium to the right ventricle and back to the lungs. Between the parts of the heart are valves, which main task is to regulate the blood flow. The motion of the heart is regulated by electrical impulses, but we will not get into these details.

2.2 The Mitral Valve

The mitral valve (MV) is the valve in the heart between the left atrium and left ventricle. It consists of two leaflets: the shorter and wider anterior, and the thinner and longer posterior leaflet. For a 2D projection of the valve, see Figure 1.

The blood flows through the MV filling the left ventricle. The filling period is called a diastole which can be divided into three parts: the early wave (E-wave), diastasis, and the atrial wave (A-wave). The E-wave is caused by active relaxation of the ventricle, causing the MV to open and blood to flow through. The A-wave is caused by contraction of the atrium and thus pushing through the remaining blood.

Between these two waves is a period of constant flow called the diastasis. When the diastole ends systole begins, which is the part of the heart cycle when the blood leaves the left ventricle through the aortic valve [7].

(6)

2.3 Numerical Hemodynamic Modeling 2 BACKGROUND

Figure 1: A 2D model of the Mitral Valve showing the anterior and posterior leaflets.

[6]

2.2.1 Mitral Valve Diseases

Mitral valve stenosis, or mitral stenosis (MS), is a condition where the orifice (open- ing) area of the MV is smaller than average. For comparison, the area is 4 to 6 cm2 for a healthy heart [8]. The condition is usually caused by rheumatic heart disease. Characteristics of MS caused by rheumatic heart disease include a fusing of the edges and thickening of the leaflets, which results in a smaller orifice. MS can also be a result of calcification or valvulitis (inflammation) [8].

Mitral regurgitation is a state where blood flows back from the left ventricle to the atrium. The condition can be caused by a rupture or restricted movements of the leaflets as for MS. MS does not usually cause regurgitation until the orifice area is below 2 cm2 [8]. No simulations of mitral valve regurgitation will be conducted in this study due to restrictions in the FEniCS-HPC HeartSolver1.

2.3 Numerical Hemodynamic Modeling

2.3.1 Navier-Stokes Equations for Moving Boundaries

Solving the Navier-Stokes equations is needed to evaluate the velocity and pressure fields of the blood flow in the left ventricle. If blood is assumed to be a Newtonian, incompressible fluid with density ρ and kinematic viscosity ν, its dynamics in a region Ω, are described by,

∂u

∂t + (u · ∇)u + 1

ρ∇p = ν∇2u, (1)

∇ · u = 0, r ∈ Ω, (2)

where u is the velocity field and p is the pressure field. Equation (1) is a momentum equation and Equation (2) is a continuity equation.

(7)

2.3 Numerical Hemodynamic Modeling 2 BACKGROUND

An Arbitrary Lagrangian-Eulerian (ALE) formulation of Navier-Stokes equations allow moving meshes [9]. This is useful in order to simulate the movement of the endocardial wall of the left ventricle. If the mesh velocity is denoted m, Navier- Stokes Equations are rewritten as (1)-(2) as

∂u

∂t + ((u − m) · ∇)u +1

ρ∇p = ν∇2u, (3)

∇ · u = 0, r ∈ Ω, (4)

u = 0, r ∈ ∂Ω. (5)

The boundary condition (5) is the no-slip condition, which says that there is no movement on the boundary due to friction and it being impenetrable. If a region allow inhomogenous Dirichlet conditions, they must also be defined. Let H ⊆ ∂Ω be a region with the Dirichlet condition f(r, t), then

u(r) = 0, r ∈ ∂Ω \ H (6)

u(r) = f(r, t), r ∈ H. (7)

If f(r, t) = 0, then the no-slip condition is everywhere and the Dirichlet conditions are homogeneous. f(r, t) 6= 0 can be interpreted as a function that describes inflow and/or outflow (depending on the sign).

2.3.2 Vorticity and Streamlines

Vorticity, ~ω is the rotation of the velocity field, defined by

~

ω = ∇ × u. (8)

Streamlines are field lines in the velocity field of a fluid and are tangent to the velocity vector at a specific time. The path a fluid particle traverses is given by the streamlines at each time. Streamlines are defined in three dimensions as

dx ux = dy

uy = dz

uz, (9)

where ux is the x-component of the velocity field at one point and time, and so on.

We remind the reader of this due to these concepts being used in the simulation results.

2.3.3 Finite Element Method

The adaptive finite element method described in [4] was used in previous studies [1]. The framework follows the scheme:

(8)

2.4 Curve Fitting 2 BACKGROUND

1. Main Navier-Stokes-solver 2. Dual Navier-Stokes-solver 3. Error estimation

4. Mesh refinement

5. Repeat until some convergence criterion is met

A tetrahedral mesh is used to represent the left ventricle of a patient specific human heart, see Figure 2. The mesh refinement can be made through identifying the longest edge of a mesh element and bisecting that edge. The computations of the mesh are parallelized (using a hybrid MPI/open mp scheme [10]), that is, different processors handle different parts of the mesh. This can lead to the bisectors not being compatible if the processors determine different bisectors for the boundary elements between them. To solve this problem, a scheme where the processors communicate with each other through voting on which edge to be bisected is introduced [10]. Note that a fixed mesh is used here.

(a) Mesh without edges, side view.

(b) Mesh with edges, side view.

(c) Mesh without edges, top view.

(d) Mesh with edges, top view.

Figure 2: Mesh used in simulations depicted during E-wave. Red: MV, Orange:

Aortic valve.

2.4 Curve Fitting

Curve fitting is useful when one is presented with discrete data and wants to trans- form it into a continuous function. Our parametric model, built on previous work [1, 2], use curve fitting to describe the motion of the short axes of the MV. This is because they are derived from measurements of the left ventricular volume which are discrete in time and only contain 27 data points.

(9)

2.4 Curve Fitting 2 BACKGROUND

2.4.1 Cubic Splines

A spline is a continuous function consisting of piece-wise polynomials, and addition- ally has a continuous first derivative. A spline can be constructed from at least two data points. A number of data points (xi, yi), 0 ≤ i ≤ N , can give rise to a spline with N piece-wise polynomials defined in each sub interval [xi, xi+1]. In the case of cubic splines, with Pi(x), xi ≤ x ≤ xi+1, it is defined as

Pi(x) = ai+ bi(x − xi) + ci(x − xi)2+ di(x − xi)3. (10) Furthermore, a natural cubic spline has the additional property that the second derivative is zero at its endpoints,

P00(x0) = P00(xN) = 0. (11) The constants ai, bi, ci and di can be found from the requirement of continuity on P (x) and Equation (11).

2.4.2 Curve Fitting with Splines

A method of curve fitting that avoid Runge’s phenomena that could emerge from polynomial interpolation is spline interpolation. However, interpolation may not yield wanted results when confronted with noisy data, and one must instead balance bias and variance. Let Yi = f (xi) + i be noisy observations where f is the true function and i are randomly picked from a distribution with zero mean and constant variance. We define the smoothing spline estimate ˆf of f that minimizes

X

i=1

(Yi− ˆf (xi))2+ λ Z

00(x)2 dx, (12)

where λ ≥ 0 is the smoothing parameter. The second term is a penalty for the roughness of the function estimate. Large λ approximates the estimate to a linear least squares estimate while small λ approximates an interpolating spline.

(10)

2.5 Simulation of a Healthy Mitral Valve 2 BACKGROUND

2.5 Simulation of a Healthy Mitral Valve

This section presents previous heart simulations of a healthy heart. This data, was supplied by Frida Svelander [1].

(a) Velocity - E-wave peak (b) Velocity - E-wave end Figure 3: A slice of the velocity field perpendicular to the MV’s long axis at two different times.

Figure 3 was made in ParaView2.The 3D-mesh was sliced parallel to the normal and perpendicular to the long axis of the MV. The figures show the velocity field in the left ventricle. The velocity magnitude is presented in m/s. Note that the area of the slice increases, a result from the ventricle increasing in volume.

2See section 3.4.3.

(11)

3 METHOD

(a) Pressure - E-wave peak (b) Pressure - E-wave end Figure 4: A slice of the pressure field perpendicular to the MV’s long axis for two different times.

Figure 4 shows the left ventricle sliced in the same way as in Figure 3 but shows the pressure field. The unit of pressure here is P a. Note that the color grading is not identical in Figure 4a and 4b, which in this case helps in visualizing the data.

3 Method

In this section, we will describe the MV model as well as introduce the modifications for modelling mitral stenosis. Lastly, we will describe the software used for this type of simulation.

3.1 Mitral Valve Design

The model used in this this study was created in previous work [1, 2]. The Mitral valve is three dimensional, but a 3D-model has not yet been implemented due to its complex geometry. Instead, a 2D-construction of the valve is used: two joined half ellipses with parallel short axes in opposite directions, placed on top of the left ventricle. Shrinking and expanding of the short axes simulate the opening and closing dynamics of the MV. The movements are derived from the heart diastole and the recorded size of the left ventricle over a diastole. The focus in this section is to describe how the model works.

(12)

3.1 Mitral Valve Design 3 METHOD

The MV model takes volume data as input. The volume data is given from echocardiac data [3] which describes how the left ventricle volume changes over time, over a diastole. This data is essential to describe when the MV should be in certain positions. For example, during the peak of the E-wave when the MV is letting blood in the most, the MV should be as open as it can be. This is not necessarily true in reality because we are not factoring in the electric impulse affecting the motion of the MV and the whole heart. In the model, the important points in time are the start, peak and end of the E-wave and the A-wave. Because the shape of the MV is a function of time, these time points must be known. By using the volume data, these timings can be extracted and then used in describing the MV movement. This procedure will be referred to as ”MV timings”.

In order to get as good MV timings as possible, the discrete volume data is interpolated using splines. The fit is done with smoothing splines, as described in section 2.4.2. However, to get the best fit, the smoothing parameter is set to the most robust value according to an algorithm based on statistical sampling of a curvature approximation. After the fit is done, the volume data can be analyzed mathematically.

In the MATLAB code, this is done in either of two ways. One is to mainly use the information about the derivative and the other information about the second derivative. The two alternatives are called ”slope” and ”curvature”. The curvature method is disregarded because it is unfinished (it does not calculate the timings for the A wave). Shown in Figure 5 are graphs of the volume data, the derivative and the second derivative denoted as ”slope” and ”curvature” respectively.

With the slope method, let v(t) be the fit of the volume data, then tsystole = t0, where dv(t0)

dt = mindv(t)

dt , 0 < t < tend, tEpeak = t2, where dv(t2)

dt = max

dv(t) dt

, tsystole < t < tend, tEstart = t1, where dv(t1)

dt = min

dv(t) dt

, tsystole < t < tEpeak, tEend = t3, where d

dt

 dv(t3) dt

d2v(t3) dt2



= 0, tEpeak < t < tend, (13) tApeak = t5, where dv(t5)

dt = maxdv(t)

dt , tEend < t < tend, tAstart = t4 6= t3, where d

dt

 dv(t4) dt

d2v(t4) dt2



= 0, tEend < t < tApeak tAend = tend, Chosen to start of systole.

However, a problem arises when there is no solution to t4 since t4 6= t3, meaning

(13)

3.1 Mitral Valve Design 3 METHOD

Figure 5: Volume data together with fit and its first and second derivative. MV timings are gathered form marked points, where the signs represent which method was used. (Star: Slope method, Circle: Curvature method)

there is only one solution to d dt

 dv(t) dt

d2v(t) dt2



= 0. (14)

This unfortunately happens with our volume data. The solution to this is to manu- ally declare where the diastasis ends. One can make a reasonable guess as to where by observing the surface area of the left ventricle over time. This declaration is presented in Figure 7, with the argument that the diastasis ends when the rate of volume change increases (A-wave). The timings needed for the movement of the MV are now found, but not the movements themselves. Geometrical data of the MV is difficult to acquire because of limited echocardiography [1]. The scaling of the leaflets’ peak position distance from the long axis relative the maximum of the leaflets’ positions are chosen with the relations shown in Table 1, for the different phases in the diastole.

Our model assumes that the leaflets are fully open at the E-wave peak (maximum open, maximum inflow) and 85 % open at the A-wave peak. Compare with Figure 6 for a visual of the anterior leaflet being in negative position at Estart. The value

−0.572 is calculated by

yEanterior

start = −yEposterior

start

Rap , (15)

where Rap is the ratio between the anterior leaflets’ max position and the posterior

(14)

3.1 Mitral Valve Design 3 METHOD

Position (ratio to max) Posterior Anterior

yEpeak 1 1

yAstart = yEend 0.7 0.5

yApeak 0.85 0.85

yAend = yEstart 0.4 -0.572

Table 1: The table shows how the positions of the posterior and anterior short axes are positioned in relation to the long axes for different time points during diastole.

The ratios are in relation to the maximum position of either short axis.

leaflets’ max position. The reason for this is that the leaflets’ starting position must be the same, since that is when the valve is closed.

Figure 6: MV model with two half ellipses with a shared long axis whose short axes move during diastole.

Figure 7: Mesh surface area versus time frame

(15)

3.1 Mitral Valve Design 3 METHOD

The timings calculated in Equation (13) and the table of positions provide a model of the dynamics of the MV. This data is also interpolated using smoothing splines3 and the same smoothing parameter that was used to fit the volume data.

Generated from the MATLAB code are two (5 × 4) matrices describing the movements of the two ellipses, one matrix for the posterior ellipsis and one matrix for the anterior. The rows describe the coefficients of a third degree polynomial (using cubic splines described in section 2.4.2) with time as the variable. We call these ramping coefficients. These polynomials scale the maximum of the ellipses’

short axes, and the splines are calculated by multiplying the matrices with the ellipses’ maximum for the short axis. The five rows each describe one part of the diastole:

1. Start of E-wave (Es) to peak of E-wave (Et) 2. Peak of E-wave to end of E-wave (Ee) 3. Diastasis

4. Start of A-wave (As) to peak of A-wave (At) 5. Peak of A-wave to end of A-wave (Ae)

Figure 8: The splines describing the movements of the ellipses’ short axes. The start and end of the splines are marked with red circles. Anterior leaflet in red and posterior leaflet in blue.

3See 2.4.2

(16)

3.1 Mitral Valve Design 3 METHOD

Figure 9: Labeling of axes in the MV model.

The time points are marked in Figure 8. The MV model is related to previous work [7], that considered the opening angels of the leaflets as parameters, instead of the short axes. That model was also based on patient specific volume data of the left ventricle. Piece-wise continuous sinusoidal functions, not cubic splines, were used to fit the data into a smooth function.

3.1.1 Model for Mitral Stenosis

Our proposed model of Mitral Stenosis (MS) is scaling the MV axes with a constant c. From the patient specific data, the following measurements were estimated from the mesh geometry:

Long axis = 14.73 mm, Short axis, posterior = 13.3 mm,

Short axis, anterior = 9.3 mm.

We relabel these as LA, SAp and SAa, according to Figure 9. Scaling these values with a constant c < 1 reduces the area of the orifice. The reduction of the area determines the severeness of the MS. The area calculation for the graphs in Figure 10 are made by:

A0 = π

2 · LA · (SAp + SAa) (16)

Introducing a scaling of each axis with a constant c, the area calculation becomes:

Ac= π

2 · cLA · (cSAp+ cSAa) ⇐⇒ Ac= π

2 · LA · (SAp+ SAa) · c2 = A0· c2 (17) Choosing c to make the area at its maximum to be smaller than 4 cm2 would put the MV orifice area below the lower limit for a healthy heart. Equation (17) yields the following relation between the axis scaling constant and the scaling of the area

c =r Ac A =√

cA with, (18)

(17)

3.2 Variation of Parameters 3 METHOD

cA = Ac

A0, (19)

where cA is the scaling constant of the original area.

3.2 Variation of Parameters

The introduction of the parameter c allows us to control the orifice area and thus the severeness of the MS. Considering that regurgitation effects supposedly starts to happen when the area < 2 cm2 [8] and the areas that are considered healthy are in the interval 4 − 6 cm2, areas we determined to study areas in the interval 2 − 3 cm2. With length of axes described earlier, reduced areas simulated correspond to the constants c in Table 2:

Areas (mm2) 300 275 250 225 200

c 0.7574 0.7252 0.6914 0.6560 0.6184

Table 2: Table over orifice areas and their corresponding scaling constants.

The constants c in Table 2 are calculated using Equation (18).

3.3 Patient Specific Data

We had access to the volume data given by echocardiography of the left ventricle in a single patient. The patient specific data gives a maximum orifice area of 522.7584 mm2 ≈ 5.23 cm2 (see Figure 10). This corresponds well to the values proposed for a healthy heart (4-6 cm2).

From the MATLAB script, described in section 3.1, using the described ”slope”

method and cubic splines, we get ramping coefficients for the MV axes over time as output. Using the provided patient specific volume data as input these ramping coefficients are:

Posterior leaflet:

−42.92 −24.21 8.49 0.4 412.66 −70.13 0 1.0

0 0 0 0.7

−1120.14 108.30 0 0.7 535.37 −141.02 0 0.85

Anterior leaflet:

−192.08 −42.69 20.89 −0.5720

687.77 −116.88 0 1.0

0 0 0 0.5

−2613.67 252.71 0 0.5

1919.89 −460.48 0 0.85

(18)

3.3 Patient Specific Data 3 METHOD

Figure 10: The orifice area over a diastole of the healthy heart from patient-specific data. The max area used for determining MS is marked by the red circle.

The coefficients are rounded for readability. For all decimals, see Appendix A. A graphic representation of the MV movements over a diastole divided into 25 time frames are shown in Figure 11.

Figure 11: A full diastole divided into 25 time points. Anterior leaflet in red and posterior in blue.

(19)

3.4 Software 3 METHOD

data:

tEs = 0.3906 s, tEt = 0.5208 s, tEe = 0.6341 s, tAs = 0.7005 s, tAt = 0.7650 s, tAe = 0.8301 s.

Splines (with rounded coefficients) to describe the motion of the posterior leaflet, with a short axis maximum of 13.3 mm:

s1p(t) = −571t3− 322t2+ 113t + 5, 0.3906 < t < 0.5208, s2p(t) = 5488t3 − 933t2+ 13, 0.5208, ≤ t < 0.6341,

s3p(t) = 9, 0.6341 ≤ t < 0.7005,

s4p(t) = −14898t3+ 1440t2+ 9, 0.7005 ≤ t < 0.7650, s5p(t) = 7120t3 − 1876t2+ 11, 0.7650 ≤ t < 0.8301

and the motion of the anterior leaflet, with a short axis maximum of 9.3 mm:

s1a(t) = −1786t3− 397t2+ 194t − 5, 0.3906 < t < 0.5208, s2a(t) = 6396t3− 1087t2+ 9, 0.5208 ≤ t < 0.6341,

s3a(t) = 5, 0.6341 ≤ t < 0.7005,

s4a(t) = −24307t3+ 2350t2+ 5, 0.7005 ≤ t < 0.7650, s5a(t) = 17855t3− 4282t2+ 8, 0.7650 ≤ t < 0.8301.

These splines are shown in Figure 8. Note that the ramping coefficients obtained are the same as in the patient specific case. An example of opening dynamics during a diastole with stenosis is shown in Figure 12, with a scaling constant of c = 0.8.

3.4 Software

3.4.1 FEniCS-HPC HeartSolver

FEniCS-HPC is a framework developed to solve PDE:s by the finite element method, built on top of DOLFIN-HPC. Applications of FEniCS include computational fluid dynamics [9], of which an example is the FEniCS-HPC HeartSolver [1, 4]. The

(20)

3.4 Software 3 METHOD

Figure 12: Example of opening with mitral valve stenosis. A c = 0.8 scaling of the ellipses’ axes to simulate a smaller orifice.

HeartSolver was developed in previous work and solves ALE formulated Navier- Stokes as described in Equation (3). We use the solver to analyze the pressure and velocity in the left ventricle. The lengths of the axes of the mitral valve as well as ramping matrices generated from MATLAB are used as input to the solver. Values for the properties of the blood was chosen the same as in previous study:

µ = 0.0027 Pa · s, ρ = 1060 kg/m3,

with ν = µ/ρ [3]. The HeartSolver was run4 on the Beskow supercomputer at PDC Center of for High Performance Computing at KTH [5].

3.4.2 DOLFIN-HPC

DOLFIN-HPC [10, 11] (high performance computing) is a branch under DOLFIN, which is an object oriented finite element library. DOLFIN-HPC is used in FEniCS- HPC HeartSolver for assembly of finite elements, handling the mesh during the par- allelized computations, and solving equations [9]. For the simulations of this study, DOLFIN-HPC was build with PETSc (for handling linear algebra), ParMETIS (for handling mesh distributions) and GTS (for intersecting geometrics).

4With the following command: srun -N 5 --ntasks-per-node=32 ./heartsolver -m data/VolumeMesh.xml -p parameters

(21)

4 RESULTS

3.4.3 DOLFIN-post and ParaView

Visualization of the simulations are produced through the open-source data visu- alization application ParaView. ParaView version 5.8.0 (Windows) is used in this paper and was released on the 18th of February 2020 [12].

DOLFIN-post is an extension of DOLFIN HPC written by Niclas Jansson. The solution files produced as output from the HeartSolver are binary and are post- processed through DOLFIN-post to vtu-files manageable by ParaView. The version of DOLFIN-post that was uploaded on the 17th of April 2016 is used here [13].

4 Results

In this section, we present the results of the simulations made with the FEniCS- HPC HeartSolver. Based on the simulations, we will make conclusions considering velocity, pressure, and vorticity.

4.1 Velocity

The velocity of the inflow during E-wave increases when the area of the orifice of the MV decreases, see Figure 13. This is due to a fix volume of blood having to flow through a smaller area. The jet of blood also reaches deeper in to the ventricle.

We see no significant difference for the different areas for diastasis, see Figure 13.

The A-wave is affected in similar ways as the E-wave: smaller orifice leads to higher inflow velocity, as shown in Figure 15.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200mm2

Figure 13: A slice of the velocity field parallel to the short axes at 540 ms (After E-wave peak). Shared color gradient.

(22)

4.2 Pressure 4 RESULTS

(a) 300 mm2 (b) 275 mm2 (c) 200 mm2 (d) 225 mm2 (e) 200 mm2

Figure 14: A slice of the velocity field parallel to the short axes at 670 ms (During diastasis). Shared color gradient.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 15: A slice of the velocity field parallel to the short axes at 780 ms (After A-wave peak). Shared color gradient.

4.2 Pressure

We see in Figure 16 the pressure distribution in the ventricle after the E-wave peak.

The area of high pressure caused by the E-wave is further down in the ventricle for ventricles with smaller MV.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 16: A slice of the pressure field parallel to the short axes at 540 ms (After E-wave peak).

(23)

4.3 Vorticity 4 RESULTS

The pressure is smaller and less concentrated at the bottom of the ventricle with smaller orifice as seen in Figure 17, but more spread out.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 17: A slice of the pressure field parallel to the short axes at 670 ms (During diastasis).

Figure 18 illustrates in all simulations a higher region of pressure just beneath the MV. This disturbs the blood from the A-wave since fluids flow from areas of high pressure to lower pressure.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 18: A slice of the pressure field parallel to the short axes at 780 ms (After A-wave peak).

4.3 Vorticity

A pressure contour is added throughout the results in this section for reference to the last section. This section illustrates the qualitative change in vorticity over time, that shows the velocity changing direction. The magnitude of the vorticity increases with decreasing area, as seen in Figure 19-21. Furthermore, the figures show all points where the magnitude of the vorticity is greater than the mean. The mean vorticity is then approximately the minimum value of the color gradient.

(24)

4.3 Vorticity 4 RESULTS

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 19: Greater than mean vorticity at 540 ms (After E-wave peak).

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 20: Greater than mean vorticity at 670 ms (During diastasis).

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2

Figure 21: Greater than mean vorticity at 780 ms (After A-wave peak).

(25)

4.4 Streamlines 4 RESULTS

4.4 Streamlines

The streamlines are additional illustrations for the reader to understand the blood flow. A pressure contour is also depicted here. The summarize this section, it is hard to draw any conclusions about the qualitative behaviour of the blood flow from these figure besides the fact that they become chaotic after the E-wave.

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2 Figure 22: Streamlines at at 540 ms (After E-wave peak).

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2 Figure 23: Streamlines at at 670 ms (During diastasis).

(26)

5 DISCUSSION

(a) 300 mm2 (b) 275 mm2 (c) 250 mm2 (d) 225 mm2 (e) 200 mm2 Figure 24: Streamlines at 780 ms (After A-wave peak).

Lastly, Figure 22 illustrate the blood flow willing to leave the ventricle through the MV, but cannot due to the Dirichlet conditions.

5 Discussion

For this final section we will summarize the results, and discuss improvements on the model based on our conclusions. We will also propose a different model for mitral stenosis.

5.1 Conclusion

As seen in chapter 4, smaller orifice caused by mitral stenosis cause greater inflow that reaches deeper in the ventricle5, greater vorticity and lower pressure throughout the ventricle. What could have clinical importance is that the decrease in overall pressure may force the heart to do more work in order to increase the volume of the ventricle. Furthermore, a lower overall pressure, visible in cases of MS, causes a lower pressure gradient between the ventricle and the aortic valve. Thus blood flowing from the left ventricle is doing so much slower, worsening blood circulation.

That is unless the heart overcompensates.

As previously mentioned, the model falls short in two regards. Firstly, the strict boundary condition for inflow, and secondly, a 2D mitral valve.

The inflow does not allow regurgitation as it is set to be a Dirichlet boundary condition: constant inflow proportional to the area of the MV. The choice of MV area was done so that regurgitation would not have to be implemented to produce realistic simulations, according to [8]. Figure 22 indicates that regurgitation might

5As seen in previous work, for example [1].

(27)

5.2 Choice of Diseases to Study 5 DISCUSSION

occur even with larger areas than 2 cm2. The approximation of the mitral valve and the simple Dirichlet boundary conditions may be causing this.

5.2 Choice of Diseases to Study

Mitral Regurgitation was not considered in this study due to being unsuitable for the given model. This disease causes blood to flow from the ventricle to the atrium be- cause the mitral valve does not close properly. The current version of FEniCS-HPC HeartSolver does not allow negative flow. Modifying the HeartSolver to account for negative flow could be a development for future studies.

5.3 Ethical Aspects

The patient-specific data was available to us through the team at KTH that con- ducted the previous study [1], which was approved to use the data through an ethical review. The patients have also given their approval of the use of their data for scientific purposes. For these types of studies the anonymization of the data is of importance, that we can not trace back what patient the data belongs to.

5.4 Further Studies

5.4.1 Shortening Axes with a Constant

Another model was considered to shorten the axes to match a MS orifice. Reducing each axis with a linear constant, k:

Long axis = LA − k mm Short axis, posterior = SAp− k mm

Short axis, anterior = SAa− k mm

Or, different k for each axes:

Long axis = LA − kLA mm Short axis, posterior = SAp− kSAp mm

Short axis, anterior = SAa− kSAa mm

It could be interesting to consider this model and compare with the scaling model used in this study. The model might be more appropriate in medical application since the constant k could represent the thickness of the calcification, for example.

(28)

REFERENCES REFERENCES

5.4.2 Negative Flow

In the simulations of this study, the flow through the MV is always in the direction from the atrium to the ventricle. One can see in Figure 22 that the blood has a tendency to flow back to the MV from the ventricle. A conclusion can be drawn that a solver that allows negative flow (from ventricle to atrium) would produce different results that could be interesting. Since the physiological MV allows for negative flow those results could be more medically accurate. This could be achieved by constructing the aortic valve and inducing a pressure gradient through the MV that does not impose Dirichlet conditions that force the blood to stay inside the left ventricle.

References

[1] Svelander et al. Towards a Framework for Patient-specific simulation of mitral valve disease. Poster Presented at: 25th Congress of the European Society of Biomechanics. Vienna. 2017.

[2] Lucor. Adaptive Cardiovascular Model Selection by Bayesian Inference. 5th In- ternational Conference on Computational and Mathematical Biomedical En- gineering. 2017.

[3] Larsson et al. Patient-Specific Left Ventricular Flow Simulations. From Transtho- racic Echocardiography:Robustness Evaluation and Validation Against Ultra- sound Doppler and Magnetic Resonance Imaging. 2017.

[4] Jeannette Hiromi Sp¨uhler. “Patient-Specific Finite Element Modeling of the Blood Flow in the Left Ventricle of a Human Heart ”. PhD thesis. KTH, 2017.

[5] PDC. Beskow. 2016. url: https : / / www . pdc . kth . se / hpc - services / computing-systems/beskow-1.737436. (accessed 27.04.2020).

[6] Wikipedia page for Mitral Valve. 2020. url: https://en.wikipedia.org/

wiki/Mitral_valve. (accessed 22.04.2020).

[7] Seo et al. “Effect of the Mitral Valve on Diastolic Flow Patterns”. In: PHYSICS OF FLUIDS 26 (2014).

[8] Serge C. Harb Brian P. Griffin. Mitral Valve Disease: a Comprehensive Review.

2017.

[9] Hoffman J. et al. “FEniCS-HPC: Coupled Multiphysics in Computational Fluid Dynamics”. In: Di Napoli E., Hermanns MA., Iliev H., Lintermann A., Peyser A. (eds) High-Performance Scientific Computing. JHPCS 2016.

Lecture Notes in Computer Science 10164 (2017).

(29)

REFERENCES REFERENCES

[10] Jansson et al. “Framework for Massively Parallel Adaptive Finite Element Computational Fluid Dynamics on Tetrahedral Meshes”. In: SIAM J. SCI.

COMPUT. 34.1 (2012), pp. C24–C41.

[11] Jansson et al. DOLFIN HPC. url: https://bitbucket.org/fenics-hpc/

dolfin-hpc/src/0.8.x/. (accessed 2.05.2020).

[12] ParaView. Official Website for ParaView. 2020. url: https://www.paraview.

org/overview/. (accessed 22.04.2020).

[13] Niclas Jansson. DOLFIN-post on BitBucket. 2016. url: https://bitbucket.

org/fenics-hpc/dolfin-post/src/master/. (accessed 22.04.2020).

(30)

A RAMPING COEFFICIENTS

Appendices

A Ramping Coefficients

The patient specific volume data as input in the MATLAB code produces the ramp- ing coefficients:

Posterior leaflet:

−42.9211851575252 −24.2098883168091 8.48802399380271 0.4

412.664219965879 −70.1250018734372 0 1.0

0 0 0 0.7

−1120.14470571677 108.304814869803 0 0.7

535.368118526369 −141.018338266496 0 0.85

Anterior leaflet:

−192.084376219415 −42.694691630016 20.88913381674 −0.5720

687.773699943132 −116.87500312239 0 1.0

0 0 0 0.5

−2613.67098000579 252.711234696208 0 0.5

1919.89199523979 −460.480618850101 0 0.85

 These were the values used as input to FEniCS-HPC Heartsolver, not the rounded

values presented in the text for readability.

References

Related documents

In a recent study, using a novel assay to differentiate serological response to HHV-6A and HHV-6B, it was demonstrated that HHV-6A serology was associated with an increased MS risk

(2000) measured a human dry skull with damping material inside and reported the resonance frequency of mechanical point impedance as 600 Hz at the posterior caudal part of

A Finite Element Model of the Human Head for Simulation of Bone

Skillnad i namn mellan attributen grapefrukt, citron och floral som alla finns med som attribut för både smak och doft borde ha varit tydligare skilda åt för att inte

Pathline visualization of the four EDV components: Direct Flow (green), Retained Inflow (yellow), Delayed Ejection Flow (Blue) and Residual Volume (Red). A-C) Healthy 56 year

Linköping University Medical Dissertations No.. Linköping University Medical

Swanson JC, Krishnamurthy G, Itoh A, Kvitting JP, Bothe W, Craig Miller D, Ingels NB, Jr. Multiple mitral leaflet contractile systems in the

The top 8 panels show, for each of the time points ❶-❽, the geometry of the anterior mitral leaflet, an outline of the mitral annulus, and the orientation of eight potential