arXiv:1704.01960v1 [physics.med-ph] 29 Mar 2017

A coupled mitral valve – left ventricle model with
fluid-structure interaction
Hao Gao
School of Mathematics and Statistics, University of Glasgow, UK

Liuyang Feng
School of Mathematics and Statistics, University of Glasgow, UK

Nan Qi
School of Mathematics and Statistics, University of Glasgow, UK

Colin Berry
Institute of Cardiovascular and Medical Science, University of Glasgow, UK

Boyce Griffith
Departments of Mathematics and Biomedical Engineering and McAllister Heart Institute,
University of North Carolina, Chapel Hill, NC, USA

Xiaoyu Luo
School of Mathematics and Statistics, University of Glasgow, UK

Abstract
Understanding the interaction between the valves and walls of the heart is important in assessing and subsequently treating heart dysfunction. With advancements in cardiac imaging, nonlinear mechanics and computational techniques, it is now possible to explore the mechanics of valve-heart interactions
using anatomically and physiologically realistic models. This study presents an
integrated model of the mitral valve (MV) coupled to the left ventricle (LV),
with the geometry derived from in vivo clinical magnetic resonance images.
Numerical simulations using this coupled MV-LV model are developed using an
immersed boundary/finite element method. The model incorporates detailed
valvular features, left ventricular contraction, nonlinear soft tissue mechanics,
and fluid-mediated interactions between the MV and LV wall. We use the model
to simulate the cardiac function from diastole to systole, and investigate how
myocardial active relaxation function affects the LV pump function. The results
of the new model agree with in vivo measurements, and demonstrate that the
diastolic filling pressure increases significantly with impaired myocardial active
relaxation to maintain the normal cardiac output. The coupled model has the

Preprint submitted to Journal of Medical Engineering & Physics

April 7, 2017

potential to advance fundamental knowledge of mechanisms underlying MV-LV
interaction, and help in risk stratification and optimization of therapies for heart
diseases.
Keywords: mitral valve, left ventricle, fluid structure interaction, immersed
boundary method, finite element method, soft tissue mechanics

2

1. Introduction
The mitral valve (MV) has a complex structure that includes two distinct
asymmetric leaflets, a mitral annulus, and chordal tendinae that connect the
leaflets to papillary muscles that attach to the wall of the left ventricle (LV).
MV dysfunction remains a major medical problem because of its close link to
cardiac dysfunctions leading to morbidity and premature mortality [1].
Computational modelling for understanding the MV mechanics promises
more effective MV repairs and replacement [2, 3, 4, 5]. Biomechanical MV
models have been developed for several decades, starting from the simplified twodimensional approximation to three-dimensional models, and to multi-physics/scale models [6, 7, 8, 9, 10, 11, 12]. Most of previous studies were based on
structural and quasi-static analysis applicable to a closed valve [13]; however,
MV function during the cardiac cycle cannot be fully assessed without modelling
the ventricular dynamics and the fluid-structure interaction (FSI) between the
MV, ventricles, and the blood flow [13, 14].
Because of the complex interactions among the MV, the sub-mitral apparatus, the heart walls, and the associated blood flow, few modelling studies have been carried out that integrate the MV and ventricles in a single
model [15, 16, 17]. Kunzelman, Einstein, and co-workers first simulated normal
and pathological mitral function [18, 19, 20] with FSI using LS-DYNA (Livermore Software Technology Corporation, Livermore, CA, USA) by putting the
MV into a straight tube. Using similar modelling approach, Lau et al. [21] compared MV dynamics with and without FSI, and they found that valvular closure
configuration is different when using the FSI MV model. Similar findings are
reported by Toma et al [22]. Over the last few years, there have also been a
number of FSI valvular models using the immersed boundary (IB) method to
study the flow across the MV [23, 24, 25]. In a series of studies, Toma [26, 22, 27]
developed a FSI MV model based on in vitro MV experimental system to study
the function of the chordal structure, and good agreement was found between
the computational model and in vitro experimental measurements. However,
none of the aforementioned MV models accounted for the MV interaction with
the LV dynamics. Indeed, Lau et al. [21] found that even with a fixed U-shaped
ventricle, the flow pattern is substantially different from that estimated using
a tubular geometry. Despite the advancements in computational modelling of
individual MV [13, 12] and LV models [28, 29, 30], it remains challenging to develop an integrated MV-LV model which includes the strong coupling between
the valvular deformation and the blood flow. Reasons for this include limited
data for model construction, difficult choices of boundary conditions, and large
computational resources required by these simulations.
Wenk et al. [15] reported a structure-only MV-LV model using LS-DYNA
that included the LV, MV, and chordae tendineae. This model was later extended to study MV stress distributions using a saddle shaped and asymmetric
mitral annuloplasty ring [16]. A more complete whole-heart model was recently
developed using a human cardiac function simulator in the Dassault Systemes’s
Living Heart project [17], which includes four ventricular chambers, cardiac
3

valves, electrophysiology, and detailed myofibre and collagen architecture. Using the same simulator, effects of different mitral annulus ring were studied by
Rausch et al. [31]. However, this simulator does not yet account for detailed
FSI.
The earliest valve-heart coupling model that includes FSI is credited to Peskin and McQueen’s pioneering work in the 1970s [32, 33, 34] using the classical
IB approach [35]. Using this same method, Yin et al. [36] investigated fluid vortices associated with the LV motion as a prescribed moving boundary. Recently,
Chandran and Kim [37] reported a prototype FSI MV dynamics in a simplified
LV chamber model during diastolic filling using an immersed interface-like approach. One of the key limitations of these coupled models is the simplified
representation of the biomechanics of the LV wall. To date, there has been no
work reported a coupled MV-LV model which has full FSI and based on realistic
geometry and experimentally-based models of soft tissue mechanics.
This study reports an integrated MV-LV model with FSI derived from in vivo
images of a healthy volunteer. Although some simplifications are made, this is
the first three-dimensional FSI MV-LV model that includes MV dynamics, LV
contraction, and experimentally constrained descriptions of nonlinear soft tissue
mechanics. This work is built on our previous models of the MV [24, 25] and
LV [38, 29]. The model is implemented using a hybrid immersed boundary
method with finite element elasticity (IB/FE) [39].
2. Methodology
2.1. IB/FE Framework
The coupled MV-LV model employs an Eulerian description for the blood,
which is modelled as a viscous incompressible fluid, along with a Lagrangian
description for the structure immersed in the fluid. The fixed physical coordinates are x = (x1 , x2 , x3 ) ∈ Ω, and the Lagrangian reference coordinate system
is X = (X1 , X2 , X3 ) ∈ U . The exterior unit normal along ∂U is N(X). Let
χ(X, t) denote the physical position of any material point X at time t, so that
χ(U, t) = Ωs (t) is the physical region occupied by the immersed structure. The
IB/FE formulation of the FSI system reads


∂u
(x, t) + u(x, t) · ∇u(x, t) = −∇p(x, t) + µ∇2 u(x, t) + f s (x, t),
(1)
ρ
∂t
∇ · u(x, t) = 0,
(2)
Z
f s (x, t) =
∇ · Ps (X, t) δ(x − χ(X, t)) dX
U
Z
−
Ps (X, t) N(X) δ(x − χ(X, t)) dA(X),
∂U

(3)
∂χ
(X, t) =
∂t

Z
u(x, t) δ(x − χ(X, t)) dx,
Ω

4

(4)

where ρ is the fluid density, µ is the fluid viscosity, u is the Eulerian velocity, p is
the Eulerian pressure, and f s is the Eulerian elastic force density. Different from
the classical IB approach [35], here the elastic force density f s is determined from
the first Piola-Kirchoff stress tensor of the immersed structure Ps as in Eq. 3.
This allows the solid deformations to be described using nonlinear soft tissue
constitutive laws. Interactions between the Lagrangian and Eulerian fields are
achieved by integral transforms with a Dirac delta function kernel δ(x) [35] in
Eqs. 3 4. For more details of the hybrid IB/FE framework, please refer to [39].
2.2. MV-LV Model Construction
A cardiac magnetic resonance (CMR) study was performed on a healthy
volunteer (male, age 28). The study was approved by the local NHS Research
Ethics Committee, and written informed consent was obtained before the CMR
scan. Twelve imaging planes along the LV outflow tract (LVOT) view were
imaged to cover the whole MV region shown in Fig. 1(a). LV geometry and
function was imaged with conventional short-axis and long-axis cine images.
The parameters for the LVOT MV cine images were: slice thickness: 3 mm with
0 gap, in-plane pixel size: 0.7×0.7 mm2 , field of view: 302 × 400 mm2 , frame
rate: 25 per cardiac cycle. Short-axis cine images covered the LV region from
the basal plane to the apex, with slice thickness: 7 mm with 3 mm gap, in-plane
pixel size: 1.3 × 1.3 mm2 , and frame rate: 25 per cardiac cycle.
The MV geometry was reconstructed from LVOT MV cine images at earlydiastole, just after the MV opens. The leaflet boundaries were manually delineated from MR images, as shown in Fig. 1(a), in which the heads of papillary
muscle and the annulus ring were identified as shown in Fig. 1(b). The MV
geometry and its sub-valvular apparatus were reconstructed using SolidWorks
(Dassault Systmes SolidWorks Corporation, Waltham, MA, USA). Because it
is difficult to see the chordal structural in the CMR, we modelled the chordae
structure using sixteen evenly distributed chordae tendineae running through
the leaflet free edges to the annulus ring, as shown in Fig. 1(c), following prior
studies [25, 24]. In a similar approach to the MV reconstruction, the LV geometry was reconstructed from the same volunteer at early-diastole by using both
the short-axis and long-axis cine images [40, 29]. Fig. 1(d) shows the inflow and
outflow tracts from one MR image. The LV wall was assembled from the short
and long axis MR images (Fig. 1(e)) to form the three dimensional reconstruction (Fig. 1(f)). The LV model was divided into four regions: the LV and the
valvular region and the inflow and the outflow tracts, as shown in Fig. 1(g).
The MV model was mounted into the inflow tract of the LV model according
to the relative positions derived from the MR images in Fig. 1(g). The left
atrium was not reconstructed but modelled as a tubular structure, the gap
between the MV annulus ring and the LV model was filled using a housing
disc structure. A three-element Windkessel model was attached to the outflow
tract of the LV model to provide physiological pressure boundary conditions
when the LV is in systolic ejection [40]. The chordae were not directly attached
to the LV wall since the papillary muscles were not modelled, similar to [25].
The myocardium has a highly layered myofibre architecture, which is usually
5

described using a fibre-sheet-normal (f , s, n) system. A rule-based method was
used to construct the myofibre orientation within the LV wall. The myofibre
angle was assumed to rotate from -60o to 60o from endocardium to epicardium,
represented by the red arrows in Fig. 1(h). In a similar way, the collagen fibres
in the MV leaflets were assumed to be circumferentially distributed, parallel
along the annulus ring, represented by the yellow arrows in Fig. 1(h).
2.3. Soft Tissue Mechanics
The total Cauchy stress (σ) in the coupled MV-LV system is
(
σ f (x, t) + σ s (x, t) for x ∈ Ωs ,
σ(x, t) =
σ f (x, t)
otherwise,

(5)

where σ f is the fluid-like stress tensor, defined as
σ f (x, t) = −pI + µ[∇u + (∇u)T ].

(6)

σ s is the solid stress tensor obtained from the nonlinear soft tissue consitutive
laws. The first Piola-Kirchhoff stress tensor Ps in Eq. 3 is related to σ s through
Ps = Jσ s F−T ,

(7)

in which F = ∂χ/∂X is the deformation gradient and J = det(F).
In the MV-LV model, we assume the structure below the LV base is contractile (Fig 1(g)), the regions above the LV basal plane, including the MV and its
apparatuses, are passive. Namely,
(
Pp + Pa below the basal plane,
s
P =
(8)
Pp
above the basal plane,
where Pa and Pp are the active and passive Piola-Kirchhoff stress tensors, respectively. The MV leaflets are modelled as an incompressible fibre-reinforced
material with the strain energy function
WMV = C1 (I1 − 3) +

av
(exp[bv (max(Ifc , 1) − 1)2 ] − 1),
2bv

(9)

in which I1 = trace(C) is the first invariant of the right Cauchy-Green deformation tensor C = FT F, Ifc = f0c · (Cf0c ) is the squared stretch along the collagen
fibre direction, and f0c denotes the collagen fibre orientation in the reference
configuration. The max() function ensures the embedded collagen network only
bears the loads when stretched, but not in compression. C1 , av , and bv are
material parameters adopted from a prior study [25] and listed in Table 1. The
passive stress tensor Pp in the MV leaflets is
Pp =

∂WMV
− C1 F−T + βs log(I3 )F−T ,
∂F
6

(10)

where I3 = det(C), and βs is the bulk modulus for ensuring the incompressibility
of immersed solid, so that the pressure-like term C1 F−T ensures the elastic stress
response is zero when F = I.
We model the chordae tendineae as the Neo-Hookean material,
Wchordae = C (I1 − 3),

(11)

where C is the shear modulus. We further assume C is much larger in systole
when the MV is closed than in diastole when the valve is opened. The much
larger value of C models the effects of papillary muscle contraction. Values of
C are listed in Table 1. Pp for the chordae tendineae is similarly derived as in
Eq. 10.
The passive response of the LV myocardium is described using the HolzapfelOgden model [41],
Wmyo =

X ai
a
exp[b(I1 − 3)] +
{exp[bi (max(I4i , 1) − 1)2 ] − 1}
2b
2bi
i=f,s

afs
{exp[bfs (I8fs )2 ] − 1}
+
2bfs

(12)

in which a, b, af , bf , as , bs , afs , bfs are the material parameters, I4f , I4s and I8fs
are the strain invariants related to the the myofibre orientations. Denoting the
myofibre direction in the reference state is f 0 and the sheet direction is s0 , we
have
I4f = f 0 · (Cf 0 ), I4s = s0 · (Cs0 ), and I8fs = f 0 · (Cs0 ).
(13)
The myocardial active stress is defined as
Pa = J T F f 0 ⊗ f 0

(14)

where T is the active tension described by the myofilament model of Niederer et
al. [42], using a set of ordinary differential equations involving the intracellular
calcium transient (Ca2+ ), sarcomere length and the active tension at the resting
sarcomere length (T req ). In our simulations, we use the same parameters as in
ref. [42], except that T req is adjusted to yield realistic contraction as the imaged
volunteer.
All the constitutive parameters in Eqs.9, 11, 12 are summarized in Table 1.
2.4. Boundary Conditions and Model Implementation
Because only the myocardium below the LV basal plane contracts, we fix the
LV basal plane along the circumferential and longitudinal displacements, but allow the radial expansion. The myocardium below the LV basal plane is left free
to move. The valvular region is assumed to be much softer than the LV region.
In diastole, a maximum displacement of 6 mm is allowed in the valvular region
using a tethering force. In systole, the valve region is gradually pulled back to
the original position. The inflow and outflow tracts are fixed. Because the MV
annulus ring are attached to a housing structure which is fixed, no additional
7

boundary conditions are applied to the MV annulus ring. Fluid boundary conditions are applied to the top planes of the inflow and outflow tracts. The function
of the aortic valve is modelled simply: the aortic valve is either fully opened
or fully closed, determined by the pressure difference between the values inside
the LV chamber and the aorta. After end-diastole, the LV region will contract
simultaneously triggered by a spatially homogeneously prescribed intracellular
Ca2+ transient [29], as shown in Fig. 3. The flow boundary conditions in a
cardiac cycle are summarized below.
• Diastolic filling: A linearly ramped pressure from 0 to a populationbased end-diastolic pressure (EDP=8 mmHg) is applied to the inflow tract
over 0.8 s, which is slightly longer than the actual diastolic duration of
the imaged volunteer (0.6 s). In diastole about 80% of diastolic filling
volume is due to the sucking effect of the left ventricle in early-diastole [43].
This negative pressure field inside the LV cavity is due to the myocardial
relaxation. We model this sucking effect using an additional pressure
loading applied to the endocardial surface, denoted as Pendo , which is
linearly ramped from 0 to 12 mmHg over 0.4 s, and then linearly decreased
to zero at end-diastole. The value of Pendo is chosen by matching the
simulated end-diastolic volume to the measured data from CMR images.
Blood flow is not allowed to move out of the LV cavity through the inflow
tract in diastole. Zero flow boundary conditions are applied to the top
plane of the outflow tract.
• Iso-volumetric contraction: Along the top plane of the inflow tract, the
EDP loading is maintained, but we allow free fluid flow in and out of the
inflow tract. Zero flow boundary conditions are retained for the outflow
tract. The duration of the iso-volumetric contraction is determined by
the myocardial contraction and ends when the aortic valve opens. The
aortic valve opens when the LV pressure is higher than the pressure in the
aorta, which is initially set to be the cuff-measured diastolic pressure in
the brachial artery, 85 mmHg.
• Systolic ejection: When the aortic valve opens, a three-element Windkessl model is coupled to the top plane of the outflow tract to provide
afterload. The volumetric flow rates across the top plane of the outflow
tract is calculated from the three-dimensional MV-LV model, and fed into
the Windkessel model [44], which returns an updated pressure for the outflow tract in the next time step. The systolic ejection phase ends when
the left ventricle cannot pump any flow through the outflow tract, and the
Windkessel model is detached.
• Iso-volumetric relaxation: Zero flow boundary conditions are applied
to both the top planes of the outflow and inflow tracts until the total cycle
ends at 1.2 seconds.
The coupled MV-LV model is immersed in a 17cm × 16cm × 16cm fluid
box. A basic time step size ∆t0 = 1.22 × 10−4 s is used in the diastolic and
8

relaxation phases, a reduced time step size (0.25 ∆t0 ) is used in the early systole
with a duration of 0.1 s, and an even smaller time step of 0.125 ∆t0 is used in
the remainder of the systolic phase. Because explicit time stepping is used in
the numerical simulations [39], we need to use a time step size small enough to
avoid numerical instabilities, particularly during the systolic phase to resolve the
highly dynamic LV deformation. The MV-LV model is implemented using the
open-source IBAMR software framework (https://github.com/IBAMR/IBAMR),
which provides an adaptive and distributed-memory parallel implementation of
the IB methods.
3. Results
Fig. 2 shows the computed volumetric flow rates across the MV and the
AV from beginning of diastole to end-systole. In diastole, the volumetric flow
rate across the MV linearly increases with Pendo , with a maximum value of
210 mL/s at 0.4 s. Diastolic filling is maintained by the increased pressure in
the inflow tract, but with decreased flow rates until end of diastole at 0.8 s. The
negative flow rate in Fig. 2 indicates the flow is entering the LV chamber. After
end-diastole, the myocardium starts to contract, and the central LV pressure
increases until it exceeds the aortic pressure (initially set to be 85 mmHg) at
0.857 s. During iso-volumetric contraction, the MV closes with a total closure
regurgitation flow of 7.2 mL, around 10% of the total filling volume, which is
comparable to the value reported by Laniado et al. [45]. There is only minor
regurgitation across the MV during systolic ejection after the iso-volumetric
contraction phase. Blood is then ejected out of the ventricle through the AV,
and the flow rate across the AV during systole reaches a peak value of 468 mL/s
(Fig. 2). The total ejection duration is 243 ms with a stroke volume of 63.2 mL.
The total blood ejected out of the LV chamber, including the regurgitation
across the MV, is 72.1 mL, which corresponds to an ejection fraction of 51%.
Fig. 3 shows the profiles of the normalized intracellular Ca2+ , LV cavity
volume, central LV pressure, and the average myocardial active tension from
diastole to systole. Until mid-diastole (0 s to 0.56 s), the central LV pressure
is negative, and the associated diastolic filling volume is around 65 mL, which
is 90% of the total diastolic filling volume. In late-diastole, the LV pressure
becomes positive. There is a delay between the myocardial active tension and
the intracellular Ca2+ profile, but the central LV pressure follows the active
tension closely throughout the cycle as shown in Fig 3.
Fig. 4 shows the deformed MV leaflets along with the corresponding CMR
cine images at early-diastole (the reference state), end-diastole, and mid-systole.
In general, the in vivo MV and LV dynamics from diastole to systole are qualitatively captured well by the coupled MV-LV model. However, a discrepancy
is observed during the diastolic filling, when the MV orifice in the model is
not opened as widely as in the CMR cine image (Fig. 4(b)). In addition, the
modelled MV leaflets have small gaps near the commissure areas even in the
fully closure state. This is partially caused by the finite size of the regularized

9

delta function at the interface and uncertainties in MV geometry reconstruction
using CMR images.
Figs. 5(a, b, c, d) show the streamlines at early-diastolic filling, late-diastolic
filling, when the MV is closing (iso-volumtric contraction), and mid-systolic
ejection when the left ventricle is ejecting. During the diastolic filling (Fig. 5(a)),
the blood flows directly through the MV into the LV chamber towards the LV
apex, in late-diastole in Fig. 5(b), the flow pattern becomes highly complex.
When iso-volmeric contraction ends, the MV is pushed back towards the left
atrium. In mid-systole, the blood is pumped out of the LV chamber through
the aortic valve into the systemic circulation, forming a strong jet as shown in
Fig. 5 (d).
The LV systolic strain related to end-diastole is shown in Fig. 6 (a), which
is negative throughout most of the region except near the basal plane, where
the LV motion is artificially constrained in the model. The average myocardial
strain along myofibre direction is -0.162±0.05. Fig. 6(b) is the fibre strain in the
MV leaflets at end-diastole, the leaflets are mostly slightly stretched during the
diastolic filling. In systole, because of the much higher pressure in the LV, the
leaflets are pushed towards the left atrium side as shown in Fig. 6(c). Near the
leaflet tip and the commissiour areas, the leaflets are highly compressed, while
in the trigons near the annulus ring, the leaflet is stretched.
From Fig. 3, one can see that the applied endocardial pressure (Pendo ) creates a negative pressure inside the LV chamber, similar to the effects of the
myocardial active relaxation. We further investigate how Pendo affects the MVLV dynamics by varying its value from 8 mmHg to 16 mmHg, and the effects
without Pendo but with an increased EDP from 8 mmHg to 20 mmHg. We observe that with an increased Pendo , the peak flow rate across the MV during the
filling phase becomes higher with more ejected volume through the aortic valve.
We also have a longer ejection duration, shorter iso-volumetric contraction time,
and higher ejection fraction as a result of increasing Pendo . On the other hand,
if we don’t apply Pendo , a much greater and nonphysiological EDP is needed for
the required ejection fraction. For example, with EDP=8 mmHg, the ejection
fraction is only 29%. Only when EDP=20 mmHg, the pump function is comparable to the case with EDP=8 mmHg and Pendo = 16 mmHg. These results are
summarized in Table 2.
4. Discussion
This study demonstrates the feasibility of integrating a MV model with a
LV model from a healthy volunteer based on in vivo CMR images. This is the
first physiologically based MV-LV model with fluid structure interaction that
includes nonlinear hyperelastic constitutive modelling of the soft tissue. The
coupled MV-LV model is used to simulate MV dynamics, LV wall deformation,
myocardial active contraction, as well as intraventricular flow. The modelling
results are in reasonable quantitative agreement with in vivo measurements
and clinical observations. For example, the peak aortic flow rate is 468 mL/s,
close to the measured peak value (498 mL/s); the ejection duration is 243 ms,
10

and the measured value is around 300 ms; the peak LV pressure is 162 mmHg,
comparable to the cuff-measured peak blood pressure 150 mmHg; the average
LV systolic strain is around -0.16, which also lies in the normal range of healthy
subjects [46].
Diastolic heart failure is usually associated with impaired myocardial relaxation and increased filling pressure [47, 48]. In this study, we model the effects of
myocardial relaxation by applying an endocardial surface pressure Pendo . Specifically we can enhance or suppress the myocardial relaxation by adjusting Pendo .
Our results in Table 2) show that, with an enhanced myocardial relaxation,
say, when Pendo ≥ 12 mmHg, there is more filling during diastole, compared to
the cases when Pendo < 12 mmHg under the same EDP. This in turn gives rise
to higher ejection fraction and stroke volume. However, if myocardial relaxation is suppressed, diastolic filling is less efficient, with subsequently smaller
ejection fraction and stroke volume. In the extreme case, when the myocardial relaxation is entirely absent, chamber volume increases by only 29.5 mL,
and ejection fraction decreases to 29%. To maintain stroke volume obtained for
Pendo =12 mmHg, EDP needs to be as high as 20 mmHg. Indeed, increased EDP
due to an impaired myocardial relaxation has been reported in a clinical study
by Zile et al. [48]. A higher EDP indicates the elevated filling pressure throughout the refilling phase. Increased filling pressure can help to maintain a normal
filling volume and ejection fraction, but runs the risks of ventricular dysfunction
in the longer term, because pump failure will occur if no other compensation
mechanism exists.
During diastole, the MV-LV model seems to yield a smaller orifice compared
to the corresponding CMR images. In our previous study [25], the MV was
mounted in a rigid straight tube, the peak diastolic filling pressure is around
10 mmHg, and the peak flow rate across the MV is comparable to the measured
value (600 mL/s). While in this coupled MV-LV model, even though with additional Pendo , the peak flow rate (200 mL/s) is much less than the measured
value. One reason is because of the extra resistance from the LV wall, which
is absent in the MV-tube model [25]. The diastolic phase can be divided into
three phasse [43]: the rapid filling, slow filling, and atrial contraction. During
rapid filling, the transvalvular flow is resulted from myocardial relaxation (the
sucking effect), which contributes to 80% of the total transvalvular flow volume.
During slow filling and atrial contraction, the left atrium needs to generate a
higher pressure to provide additional filling. In the coupled MV-LV model,
the ramped pressure in the top plane of the inflow tract during late-diastole is
related to the atrial contraction, and during this time, only 10% of the total
transvalvular flow occurs. However, the peak flow rate in rapid filling phase
is much lower compared to the measured value, which suggests the myocardial
relaxation would be much stronger.
In a series of studies based on in vitro µCT experiments, Toma [26, 27, 22]
suggested that MV models with simplified chordal structure would not compare
well with experimental data, and that a subject-specific 3D chordal structure
is necessary. This may explain some of the discrepancies we observed here. A
simplified chordal structure is used in this study because we are unable to re11

construct the chordal structure from the CMR data. CT imaging may allow the
chordae reconstruction but it comes with radiation risk. Patient-specific chordal
structure in the coupled MV-LV model would require further improvements of
in vivo imaging techniques.
Several other limitations in the model may also contribute to the discrepancies. These include the uncertainty of patient-specific parameter identification,
uncertainties in MV geometry reconstruction from CMR images, the passive
response assumption around the annulus ring and the valvular region of the LV
model, and the lack of pre-strain effects. Studies addressing these issues are
already under way. We expect that further improvement in personalized modelling and more efficient high performance computing would make the modelling
more physiologically detailed yet fast enough for applications in risk stratification and optimization of therapies in heart diseases.
5. Conclusion
We have developed a first fully coupled MV-LV model that includes fluidstructure interaction as well as experimentally constrained descriptions of the
soft tissue mechanics. The model geometry is derived from in vivo magnetic resonance images of a healthy volunteer. It incorporates three-dimensional finite
element representations of the MV leaflets, sub-valvular apparatus, and the LV
geometry. Fibre-reinforced hyperelastic constitutive laws are used to describe
the passive response of the soft tissues, and the myocardial active contraction is
also modelled. The developed MV-LV model is used to simulate MV dynamics,
LV wall deformation, and ventricular flow throughout the cardiac cycle. Despite
several modelling limitations, most of the results agree with in vivo measurements. We find that with impaired myocardial active relaxation, the diastolic
filling pressure needs to increase significantly in order to maintain a normal cardiac output, consistent with clinical observations. The model thereby represents
a further step towards a whole-heart multiphysics modelling with a target for
clinical applications.
Acknowledgement
We are grateful for the funding from the UK EPSRC (EP/N014642/1, and
EP/I029990/1) and the British Heart Foundation (PG/14/64/31043), and the
National Natural Science Foundation of China (No. 11471261). In addition,
Feng received the China Scholarship Council Studentship and the Fee Waiver
Programme at the University of Glasgow, Luo is funded by a Leverhulme Trust
Fellowship (RF-2015-510), and Griffith is supported by the National Science
Foundation (NSF award ACI 1450327) and the National Institutes of Health
(NIH award HL117063)
Conflict interests
The authors have no conflicts of interest.
12

References
References
[1] A. S. Go, D. Mozaffarian, V. L. Roger, E. J. Benjamin, J. D. Berry, M. J.
Blaha, S. Dai, E. S. Ford, C. S. Fox, S. Franco, et al., Heart disease and
stroke statistics-2014 update, Circulation 129 (3).
[2] M. S. Sacks, W. David Merryman, D. E. Schmidt, On the biomechanics of
heart valve function, Journal of biomechanics 42 (12) (2009) 1804–1824.
[3] E. Votta, T. B. Le, M. Stevanella, L. Fusini, E. G. Caiani, A. Redaelli,
F. Sotiropoulos, Toward patient-specific simulations of cardiac valves:
State-of-the-art and future directions, Journal of biomechanics 46 (2)
(2013) 217–228.
[4] W. Sun, C. Martin, T. Pham, Computational modeling of cardiac valve
function and intervention, Annual review of biomedical engineering 16
(2014) 53–76.
[5] A. Kheradvar, E. M. Groves, L. P. Dasi, S. H. Alavi, R. Tranquillo, K. J.
Grande-Allen, C. A. Simmons, B. Griffith, A. Falahatpisheh, C. J. Goergen,
et al., Emerging trends in heart valve engineering: Part i. solutions for
future, Annals of biomedical engineering 43 (4) (2015) 833–843.
[6] K. S. Kunzelman, R. Cochran, Stress/strain characteristics of porcine mitral valve tissue: parallel versus perpendicular collagen orientation, Journal
of cardiac surgery 7 (1) (1992) 71–78.
[7] S. K. Dahl, J. Vierendeels, J. Degroote, S. Annerel, L. R. Hellevik,
B. Skallerud, Fsi simulation of asymmetric mitral valve dynamics during
diastolic filling, Computer methods in biomechanics and biomedical engineering 15 (2) (2012) 121–130.
[8] E. J. Weinberg, D. Shahmirzadi, M. R. Kaazempur Mofrad, On the multiscale modeling of heart valve biomechanics in health and disease, Biomechanics and modeling in mechanobiology 9 (4) (2010) 373–387.
[9] Q. Wang, W. Sun, Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans,
Annals of biomedical engineering 41 (1) (2013) 142–153.
[10] V. Prot, B. Skallerud, Nonlinear solid finite element analysis of mitral valves
with heterogeneous leaflet layers, Computational Mechanics 43 (3) (2009)
353–368.
[11] M. Stevanella, F. Maffessanti, C. A. Conti, E. Votta, A. Arnoldi, M. Lombardi, O. Parodi, E. G. Caiani, A. Redaelli, Mitral valve patient-specific
finite element modeling from cardiac mri: Application to an annuloplasty
procedure, Cardiovascular Engineering and Technology 2 (2) (2011) 66–76.
13

[12] C.-H. Lee, C. A. Carruthers, S. Ayoub, R. C. Gorman, J. H. Gorman, M. S.
Sacks, Quantification and simulation of layer-specific mitral valve interstitial cells deformation under physiological loading, Journal of theoretical
biology 373 (2015) 26–39.
[13] D. R. Einstein, F. Del Pin, X. Jiao, A. P. Kuprat, J. P. Carson, K. S. Kunzelman, R. P. Cochran, J. M. Guccione, M. B. Ratcliffe, Fluid–structure
interactions of the mitral valve and left heart: comprehensive strategies,
past, present and future, International Journal for Numerical Methods in
Biomedical Engineering 26 (3-4) (2010) 348–380.
[14] H. Gao, N. Qi, L. Feng, X. Ma, M. Danton, C. Berry, X. Luo, Modelling
mitral valvular dynamics–current trend and future directions, International
Journal for Numerical Methods in Biomedical Engineeringdoi:10.1002/
cnm.2858.
[15] J. F. Wenk, Z. Zhang, G. Cheng, D. Malhotra, G. Acevedo-Bolton,
M. Burger, T. Suzuki, D. A. Saloner, A. W. Wallace, J. M. Guccione,
et al., First finite element model of the left ventricle with mitral valve:
insights into ischemic mitral regurgitation, The Annals of thoracic surgery
89 (5) (2010) 1546–1553.
[16] V. M. Wong, J. F. Wenk, Z. Zhang, G. Cheng, G. Acevedo-Bolton,
M. Burger, D. A. Saloner, A. W. Wallace, J. M. Guccione, M. B. Ratcliffe, et al., The effect of mitral annuloplasty shape in ischemic mitral
regurgitation: a finite element simulation, The Annals of thoracic surgery
93 (3) (2012) 776–782.
[17] B. Baillargeon, I. Costa, J. R. Leach, L. C. Lee, M. Genet, A. Toutain, J. F.
Wenk, M. K. Rausch, N. Rebelo, G. Acevedo-Bolton, et al., Human cardiac function simulator for the optimal design of a novel annuloplasty ring
with a sub-valvular element for correction of ischemic mitral regurgitation,
Cardiovascular engineering and technology 6 (2) (2015) 105–116.
[18] D. R. Einstein, P. Reinhall, M. Nicosia, R. P. Cochran, K. Kunzelman, Dynamic finite element implementation of nonlinear, anisotropic hyperelastic
biological membranes, Computer Methods in Biomechanics and Biomedical
Engineering 6 (1) (2003) 33–44.
[19] D. R. Einstein, K. S. Kunzelman, P. G. Reinhall, M. A. Nicosia, R. P.
Cochran, Non-linear fluid-coupled computational model of the mitral valve,
Journal of Heart Valve Disease 14 (3) (2005) 376–385.
[20] K. Kunzelman, D. R. Einstein, R. Cochran, Fluid–structure interaction
models of the mitral valve: function in normal and pathological states,
Philosophical Transactions of the Royal Society B: Biological Sciences
362 (1484) (2007) 1393–1406.

14

[21] K. Lau, V. Diaz, P. Scambler, G. Burriesci, Mitral valve dynamics in structural and fluid–structure interaction models, Medical engineering & physics
32 (9) (2010) 1057–1064.
[22] M. Toma, D. R. Einstein, C. H. Bloodworth, R. P. Cochran, A. P. Yoganathan, K. S. Kunzelman, Fluid–structure interaction and structural
analyses using a comprehensive mitral valve model with 3d chordal structure, International Journal for Numerical Methods in Biomedical Engineeringdoi:10.1002/cnm.2815.
[23] P. N. Watton, X. Y. Luo, M. Yin, G. M. Bernacca, D. J. Wheatley, Effect
of ventricle motion on the dynamic behaviour of chorded mitral valves,
Journal of Fluids and Structures 24 (1) (2008) 58–74.
[24] X. Ma, H. Gao, B. E. Griffith, C. Berry, X. Luo, Image-based fluid–
structure interaction model of the human mitral valve, Computers & Fluids
71 (2013) 417–425.
[25] H. Gao, N. Ma, X.and Qi, C. Berry, B. E. Griffith, X. Y. Luo, A finite strain
nonlinear human mitral valve model with fluid-structure interaction, International journal for numerical methods in biomedical engineering 30 (12)
(2014) 1597–1613.
[26] M. Toma, M. Ø. Jensen, D. R. Einstein, A. P. Yoganathan, R. P. Cochran,
K. S. Kunzelman, Fluid–structure interaction analysis of papillary muscle
forces using a comprehensive mitral valve model with 3d chordal structure,
Annals of biomedical engineering 44 (4) (2016) 942–953.
[27] M. Toma, C. H. Bloodworth, E. L. Pierce, D. R. Einstein, R. P. Cochran,
A. P. Yoganathan, K. S. Kunzelman, Fluid-structure interaction analysis
of ruptured mitral chordae tendineae, Annals of Biomedical Engineering
(2016) 1–13doi:10.1007/s10439-016-1727-y.
[28] M. P. Nash, P. J. Hunter, Computational mechanics of the heart, Journal
of elasticity and the physical science of solids 61 (1-3) (2000) 113–141.
[29] W. W. Chen, H. Gao, X. Y. Luo, N. A. Hill, Study of cardiovascular function using a coupled left ventricle and systemic circulation model, Journal
of Biomechanics 49 (12) (2016) 2445–2454. doi:10.1016/j.jbiomech.
2016.03.009.
[30] A. Quarteroni, T. Lassila, S. Rossi, R. Ruiz-Baier, Integrated heart–
coupling multiscale and multiphysics models for the simulation of the cardiac function, Computer Methods in Applied Mechanics and Engineering
314 (2016) 345–407.
[31] M. K. Rausch, A. M. Zöllner, M. Genet, B. Baillargeon, W. Bothe, E. Kuhl,
A virtual sizing tool for mitral valve annuloplasty, International Journal for
Numerical Methods in Biomedical Engineeringdoi:10.1002/cnm.2788.

15

[32] C. S. Peskin, Flow patterns around heart valves: a numerical method,
Journal of computational physics 10 (2) (1972) 252–271.
[33] D. M. McQueen, C. S. Peskin, E. L. Yellin, Fluid dynamics of the mitral
valve: physiological aspects of a mathematical model, American Journal of
Physiology-Heart and Circulatory Physiology 242 (6) (1982) H1095–H1110.
[34] C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of
computational physics 25 (3) (1977) 220–252.
[35] C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002)
479–517.
[36] M. Yin, X. Y. Luo, T. J. Wang, P. N. Watton, Effects of flow vortex
on a chorded mitral valve in the left ventricle, International Journal for
Numerical Methods in Biomedical Engineering 26 (3-4) (2010) 381–404.
[37] K. B. Chandran, H. Kim, Computational mitral valve evaluation and potential clinical applications, Annals of biomedical engineering 43 (6) (2015)
1348–1362.
[38] H. Gao, H. Wang, C. Berry, X. Y. Luo, B. E. Griffith, Quasi-static imagebased immersed boundary-finite element model of left ventricle under diastolic loading, International journal for numerical methods in biomedical
engineeringdoi:10.1002/cnm.2652.
[39] B. E. Griffith, X. Y. Luo, Hybrid finite difference/finite element
version of the immersed boundary method, eprint from arXiv (url:
https://arxiv.org/abs/1612.05916).
[40] H. Gao, C. Berry, X. Y. Luo, Image-derived human left ventricular modelling with fluid-structure interaction, in: Functional Imaging and Modeling
of the Heart, Springer, 2015, pp. 321–329.
[41] G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization,
Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 367 (1902) (2009) 3445–3475.
[42] S. Niederer, P. Hunter, N. Smith, A quantitative analysis of cardiac myocyte relaxation: a simulation study, Biophysical journal 90 (5) (2006)
1697–1722.
[43] R. A. Nishimura, A. J. Tajik, Evaluation of diastolic filling of left ventricle
in health and disease: Doppler echocardiography is the clinician?s rosetta
stone, Journal of the American College of Cardiology 30 (1) (1997) 8–18.
[44] B. E. Griffith, Immersed boundary model of aortic heart valve dynamics
with physiological driving and loading conditions, International Journal for
Numerical Methods in Biomedical Engineering 28 (3) (2012) 317–345.
16

[45] S. Laniado, E. Yellin, M. Kotler, L. Levy, J. Stadler, R. Terdiman, A study
of the dynamic relations between the mitral valve echogram and phasic
mitral flow, Circulation 51 (1) (1975) 104–113.
[46] K. Mangion, H. Gao, C. McComb, D. Carrick, G. Clerfond, X. Zhong,
X. Luo, C. E. Haig, C. Berry, A novel method for estimating myocardial
strain: Assessment of deformation tracking against reference magnetic resonance methods in healthy volunteers, Scientific Reports 6 (2016) 38774.
doi:10.1038/srep38774.
[47] I. Hay, J. Rich, P. Ferber, D. Burkhoff, M. S. Maurer, Role of impaired
myocardial relaxation in the production of elevated left ventricular filling
pressure, American Journal of Physiology-Heart and Circulatory Physiology 288 (3) (2005) H1203–H1208.
[48] M. R. Zile, C. F. Baicu, W. H. Gaasch, Diastolic heart failureabnormalities
in active relaxation and passive stiffness of the left ventricle, New England
Journal of Medicine 350 (19) (2004) 1953–1959.

17

MV leaflets
C1 (kPa)
Anterior
17.4
Posterior
10.2
Chordae
C (kPa)
systole
9000
diastole
540
Myocardium
passive
a (kPa)
b
af (kPa)
0.24
5.08
1.46
active
T req = 225 kPa

av (kPa)
31.3
50.0

bf
4.15

as (kPa)
0.87

bv (kPa)
55.93
63.48

bs
1.6

a8fs (kPa)
0.3

Table 1: Material parameter values for MV leaftlets, chordae and the myocardium

18

b8fs
1.3

Cases (mmHg)
EDP=8, Pendo =8
EDP=8, Pendo =10
EDP=8, Pendo =12
EDP=8, Pendo =14
EDP=8, Pendo =16
EDP=8, Pendo =0
EDP=12,Pendo =0
EDP=14,Pendo =0
EDP=16,Pendo =0
EDP=18,Pendo =0
EDP=20,Pendo =0

tiso-con (ms)
60
58
57
55
54
75
64
61
58
56
55

tejection (ms)
227
237
243
251
256
174
213
226
243
251
262

ejection
VLV
(mL)
52
57.6
63.2
67.8
72.3
20.8
41.0
50.6
61.8
68.5
75.7

filling
VMV
(mL)
60.6
65.9
72.1
76.8
81.3
29.5
50.9
59.8
71.9
79.3
86.2

peak
FMV
(mL/s)
412.93
442.60
468.41
486.84
503.76
209.54
343.47
406.81
459.64
486.58
511.16

Table 2: Effects of EDP and the endocardial pressure loading (Pendo ) on MV and LV dynamics.

19

LVEF(%)
47%
49%
51%
53%
54%
29%
42%
47%
51%
54%
55%

(a)

(b)

(d)

(g)

(c)

(e)

(f)

(h)

Figure 1: The CMR-derived MV-LV model. (a) The MV leaflets were segmented from a stack
of MR images of a volunteer at early-diastole, (b) positions of the papillary muscle heads and
the annulus ring, (c) reconstructed MV geometry with chordae, (d) a MR image showing the
LV and location of the outflow tract (AV) and inflow tract (MV), (e) the LV wall delineation
from short and long axis MR images, (f) the reconstructed LV model, in which the LV model
is divided into four part: the LV region bellow the LV base, the valvular region, and the inflow
and outflow tracts, (g) the rule-based fibre orientations in the LV and the MV, and (h) the
coupled MV-LV model.

20

Figure 2: Flow rates across the AV and MV from diastole to systole. Diastolic phase: 0 s to
0.8 s; Systolic phase: 0.8 s and onwards. Positive flow rate means the blood flows out of the
LV chamber.

21

Figure 3: Normalized intracellular Ca2 +, LV cavity volume, central LV pressure and average
myocardial active tension. All curves are normalized to their own maximum values, which are:
1µMol for Ca2 +, 145 mL for LV cavity volume, 162 mmHg for central LV pressure, 96.3 kPa
for average myocardial active tension.

22

(a)

(b)

(c)
Figure 4: Comparisons between the MV and LV structures at (a) reference configuration, (b)
end-diastole, and (c) end-systole, and the corresponding CMR cine images (left). Coloured
by the displacement magnitude.

23

(a)

(b)

(c)

(d)

Figure 5: Streamlines in the MV-LV model at early-diastolic filling (a), late-diastolic filling
(b), when isovolumtric contraction ends (c), and at the mid-systole. Streamline are colored
by velocity magnitude, the LV wall and MV are colored by the displacement magnitude. Red
: high; blue: low

24

(a)

(b)

(c)

Figure 6: Distributions of fibre strain in the left ventricle at end-systole (a), in the MV at
end-diastole (b) and end-systole (c).

25

