A coupled mitral valve—left ventricle model with fluid–structure interaction

Understanding the interaction between the valves and walls of the heart is important in assessing and subsequently treating heart dysfunction. 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 cardiac function from diastole to systole. Numerically predicted LV pump function agrees well with in vivo data of the imaged healthy volunteer, including the peak aortic flow rate, the systolic ejection duration, and the LV ejection fraction. In vivo MV dynamics are qualitatively captured. We further demonstrate that the diastolic filling pressure increases significantly with impaired myocardial active relaxation to maintain a normal cardiac output. This is consistent with clinical observations. The coupled model has the potential to advance our fundamental knowledge of mechanisms underlying MV-LV interaction, and help in risk stratification and optimisation of therapies for heart diseases.


Introduction
The mitral valve (MV) has a complex structure that includes two distinct asymmetric leaflets, a mitral annulus, and chordae tendineae 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 dysfunction leading to morbidity and premature mortality [1].
Computational modelling for understanding 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 two-dimensional approximation to three-dimensional models, and to multi-physics/-scale models [6][7][8][9][10][11][12]. Most 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 between the MV, the sub-mitral apparatus, the heart walls, and the associated blood flow, very few modelling studies have been carried out that integrate the MV and ventricles in a single model [15][16][17]. Kunzelman, Einstein, and coworkers first simulated normal and pathological mitral function [18][19][20] with FSI using LS-DYNA (Livermore Software Technology Corporation, Livermore, CA, USA) by putting a MV into a straight tube. Using a 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 were reported by Toma et al. [22]. Over the last few years, there have also been a number of FSI valvular models using immersed boundary (IB) method to study the flow across the MV [23][24][25]. In a series of studies, Toma et al. [22,26,27] developed an 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 was substantially different from that estimated using a tubular geometry. Despite the advancements in computational modelling of individual MV [12,13] and LV models [28][29][30], it remains challenging to develop an integrated MV-LV model that 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. motion as a prescribed moving boundary. Recently, Chan-dran 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 FSI model 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 [29,38]. The model is implemented using a hybrid immersed boundary method with finite element elasticity (IB/FE) [39].

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 were 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 mm 2 ; field of view: 302 × 400mm 2 ; 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 mm 2 ; frame rate: 25 per cardiac cycle.
The MV geometry was reconstructed from LVOT MV cine images at early-diastole, 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 Systèmes 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 [24,25]. 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 [29,40]. Fig. 1(d) shows the inflow and outflow tracts in one MR image. The LV wall was assembled from the short-/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 wall, 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, and 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, as in our previous study [25]. The myocardium has a highly layered myofibre architecture, which is usually 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 −60° to 60° 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).

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. The fixed physical coordinates are x = (x 1 ,x 2 ,x 3 ) ∈ Ω, and the Lagrangian reference coordinate system is X = (X 1 ,X 2 ,X 3 ) ∈ 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 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, which is determined from the deformed immersed structure and its active contraction. Note the IB formulation employs a common momentum equation for both the fluid and the solid, in which the additional solid stresses are accounted for by f s in Eq. (1). Interactions between the Lagrangian and Eulerian fields are achieved by integral transforms with a Dirac delta function kernel δ(x) [35], in Eqs.
(3) and (4). Different from the classical IB approach [35], here the elastic force density f s is to be described using a nonlinear hyperelastic constitutive law; see Section 2.3 below. For more details of the hybrid IB/FE framework, readers are referred to the paper by Griffith and Luo [39].

Soft tissue mechanics
The total Cauchy stress (σ) in the coupled MV-LV system is where σ f is the fluid-like stress tensor, defined as σ s is the solid stress tensor obtained from the nonlinear soft tissue constitutive laws. The first in which = ∂χ/ ∂X is the deformation gradient and J = det( ).
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, ℙ s = ℙ p + ℙ a below the basal plane, ℙ p above the basal plane, where ℙ a and ℙ p 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 in which I 1 = trace(ℂ) is the first invariant of the right Cauchy-Green deformation tensor c is the squared stretch along the collagen fibre direction, and f 0 c denotes the collagen fibre orientation in the reference configuration. The max (•) function ensures the embedded collagen network only bears loads when stretched, but not in compression. C 1 , a v , and b v are material parameters adopted from a prior study [25] and listed in Table 1. The passive stress tensor ℙ p in the MV leaflets is then where I 3 = det(ℂ), and β s is the bulk modulus for ensuring the incompressibility of immersed solid, and the pressure-like term C 1 −T ensures the elastic stress response is zero when = .
We model the chordae tendineae as a neo-Hookean material, 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 to account for the effects of papillary muscle contraction. Values of C are listed in Table 1. ℙ p for the chordae tendineae is similarly derived as in Eq. (10).
The passive response of the LV myocardium is described using the Holzapfel-Ogden model [41], in which a, b, a f , b f , a s , b s , a fS , b fs are the material parameters, I 4f, I 4s and I 8fs are the strain invariants related to the myofibre orientations. Denoting the myofibre direction in the reference state f 0 and the sheet direction by s 0 we have The myocardial active stress is defined as 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 (Ca 2+ ), sarcomere length, and the active tension at the resting sarcomere length (T req ).
The constitutive parameters in Eqs. (9), (11) and (12), summarised in Table 1, are obtained from our previous studies [25,29] which showed that the simulated MV and LV dynamics matched the in vivo measurements well. In the active contraction model (Eq. (14)), we adopted the same parameters as used by Niederer et al. [42], except that T req is set to 225 kPa to match the contractions observed in the imaged volunteer.

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 is attached to a housing structure which is fixed, no additional 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 simply modelled as: 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 enddiastole, the LV region will contract simultaneously triggered by a spatially homogeneously prescribed intracellular Ca 2+ transient [29], as shown in Fig.3. The flow boundary conditions in a cardiac cycle are summarised below.
• Diastolic filling: A linearly ramped pressure from 0 to a population-based enddiastolic 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 P endo , 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 P endo 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, and the flow is controlled by the MV leaflet dynamics. Note during the iso-volumetric contraction, regurgitation may occur due to the MV closure action and the dysfunction of the MV apparatus. However, a small backflow before the MV is fully closed may be deemed normal [44]. Zero flow boundary conditions are retained for the outflow tract. The duration of the isovolumetric 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 cuffmeasured 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 [45], 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.  Fig. 2 shows the computed volumetric flow rates through the MV and the AV from beginning of diastole to end-systole. In diastole, the volumetric flow rate through the MV linearly increases with increased pressure loading in the endocardial surface, with a maximum value of 210mL/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. [44]. 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 through the AV during systole reaches a peak value of 468 mL/s (Fig. 2) (the CMR measured value is 498 mL/s). The total ejection duration is 243 ms (the measured duration is 300 ms) with a stroke volume of 63.2 mL. The total blood ejected out of the LV chamber, including the regurgitation through the MV, is 72.1 mL. This corresponds to an ejection fraction of 51%, and is close to the CMR measured value of 57%. Fig. 3 shows the profiles of the normalised intracellular Ca 2+ , LV cavity volume, central LV pressure, and the average myocardial active tension from diastole to systole. Until middiastole (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 Ca 2 + profile, but the central LV pressure follows the active tension closely throughout the cycle as shown in Fig. 3. The peak systolic LV pressure is 162 mmHg, comparable to the cuff-measured peak blood pressure 150 mmHg. Fig. 4(a-d) shows the streamlines at early-diastolic filling, late-diastolic filling, when the MV is closing (iso-volumetric contraction), and mid-systolic ejection when the left ventricle is ejecting, During the diastolic filling ( Fig. 4(a)), the blood flows directly through the MV into the LV chamber towards the LV apex, in late-diastole in Fig. 4(b), the flow pattern becomes highly complex. When iso-volumetric 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. 4(d).

MV and myocardial dynamics
Fig . 5 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. 5(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 regularised delta function at the interface and uncertainties in MV geometry reconstruction using CMR images.
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, which lies in the normal range of healthy subjects [46]. 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.

Effects of P endo on pump function
From Fig. 3, it is clear that that the applied endocardial pressure (P endo ) creates a negative pressure inside the LV chamber, similar to the effects of the myocardial active relaxation. We further investigate how P endo affects the MV-LV dynamics by varying its value from 8 mmHg to 16 mmHg, and the effects without P endo but with an increased EDP from 8 mmHg to 20 mmHg. We observe that with an increased P endo , 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 P endo . On the other hand, if we do not apply P endo, a much greater and non-physiological 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 P endo = 16 mmHg. These results are summarised in Table 2.

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, such as the peak aortic flow rate (468mL/s vs. 498mL/s), the ejection duration (243 ms vs. 300 ms), peak cuff-measured systolic pressure (162 mmHg vs. 150 mmHg), and LV ejection fraction (51% vs. 57%). Note that the above comparisons are made with computational vs. imaged values.
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 P endo. Specifically, we can enhance or suppress the myocardial relaxation by adjusting P endo . Our results in Table 2 show that, with an enhanced myocardial relaxation, say, when P endo ≥ 12 mmHg, there is more filling during diastole, compared to the cases when P endo < 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 P endo =12 mmHg, EDP needs to be as high as 20 mmHg. Indeed, increased EDP resulted from 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 P endo, the peak flow rate (200mL/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 phases [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 that the myocardial relaxation would be much stronger.
In a series of studies based on in vitro μCT experiments, Toma et al. [22,26,27] 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 observe here. A simplified chordal structure is used in this study because we are unable to reconstruct the chordal structure from the CMR data. CT imaging may allow the chordae reconstruction but it comes with radiation risk. Patientspecific 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, the uncertainty 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 personalised modelling and more efficient high performance computing would make the modelling more physiologically detailed yet fast enough for applications in risk stratification and optimisation of therapies in heart diseases.

Conclusion
Interaction between the mitral valve and the ventricular wall plays an essential role in cardiac pump function. In this study, we have developed a first fully coupled MV-LV model that includes fluid-structure interaction as well as soft tissue mechanics with model parameters determined from in vivo image data. The model geometry is derived from in vivo magnetic resonance images of a healthy volunteer, and 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 a myofilament model is used to model the myocardial active contraction. Our results show that the developed MV-LV model can simulate MV-LV interaction with good agreements, including the peak aortic flow rate, the peak systemic pressure, the systolic ejection duration and the LV ejection fraction, with in vivo measurements despite several modelling limitations. We further find that with impaired myocardial active relaxation, the diastolic filling pressure needs to increase significantly in order to maintain a normal cardiac output, which is consistent with clinical observations in patients with impaired myocardial relaxation. This model thereby represents a step towards a whole-heart multiphysics modelling with a target for clinical applications.     Table 1 Material parameter values for MV leaflets, chordae and the myocardium.