Some Effects of Different Constitutive Laws on FSI Simulation for the Mitral Valve

In this paper, three different constitutive laws for mitral leaflets and two laws for chordae tendineae are selected to study their effects on mitral valve dynamics with fluid-structure interaction. We first fit these three mitral leaflet constitutive laws and two chordae tendineae laws with experimental data. The fluid-structure interaction is implemented in an immersed boundary framework with finite element extension for solid, that is the hybrid immersed boundary/finite element(IB/FE) method. We specifically compare the fluid-structure results of different constitutive laws since fluid-structure interaction is the physiological loading environment. This allows us to look at the peak jet velocity, the closure regurgitation volume, and the orifice area. Our numerical results show that different constitutive laws can affect mitral valve dynamics, such as the transvalvular flow rate, closure regurgitation and the orifice area, while the differences in fiber strain and stress are insignificant because all leaflet constitutive laws are fitted to the same set of experimental data. In addition, when an exponential constitutive law of chordae tendineae is used, a lower closure regurgitation flow is observed compared to that of a linear material model. In conclusion, combining numerical dynamic simulations and static experimental tests, we are able to identify suitable constitutive laws for dynamic behaviour of mitral leaflets and chordae under physiological conditions.


Results
The effects of different constitutive laws of MV leaflets. With the same linear elastic constitutive law of chordae tendineae, cases M1, M2 and M3 are used to describe the mechanical properties of valve leaflets. For all three cases, the MV is fully opened at t = 0.1 s with a pressure gradient of 10 mmHg; at t = 0.22 s, the mitral valve is just closed with a pressure difference of around 80 mmHg; and at t = 0.35 s, the MV is fully loaded, the pressure gradient reaches 150 mmHg. Table 1 shows the average and maximum displacements with three constitutive laws at fully-opened and fully-loaded states. We can see that there are some differences among the three cases. For example, the average displacement from case M3 is the largest when fully opened, while that of case M1 is the largest when fully loaded. This indicates that different leaflets constitutive laws would affect MV dynamics with FSI. Figure 1 shows the velocity fields at fully-opened and just-closed states using the three constitutive laws of leaflets, respectively. It can be found that at fully-opened state ( Fig. 1(a,c,e)), a strong jet is formed toward the outlet (the left ventricle side). When the MV is just closed ( Fig. 1(b,d,f)), the MV leaflets prevent further blood flowing back into the left atrium side with a clear closure regurgitation, especially in cases M1 and M3. Although the general flow fields in the three cases are similar, there are some minor differences. For example, the jet in case M1 at fully-opened state is stronger than the other two cases, while the closure regurgitation flow is stronger in case M3 than other two cases. The peak jet velocity at different times are given in Table 2. Slightly lower peak velocities can be found in case M2. Figures 2 and 3 show the fiber strain and stress distributions. When the MV leaflets are fully opened, both the anterior and posterior leaflets are stretched along the fiber direction with low stress due to low diastolic loading. With increased transvalvular pressure, the valve leaflets start to close. When the valve is fully loaded, regions near commissures are highly compressed in cases M1 and M3, but not in case M2. The MV leaflets are further pushed towards the left atrium side with increased pressure, better closer configurations are achieved in cases M1 and M3 with smaller orifice area compared to case M2 at fully-loaded state. When the MV is closed, all cases show high stress concentration near the annulus region, the stress level in case M2 is higher than the other two cases. In order to further analyze fiber strain and stress of leaflets, we define three local regions at the anterior leaflet: two trigones and one belly region, as shown in Fig. 4. The average fiber stress and strain of these regions are summarized in Tables 3 and 4, respectively. When the MV is fully opened, the stress in the belly region is greater than that in the trigone regions; when the MV is just closed, the stress of leaflets begin to increase with much higher  values in the trigones than that in the belly region; the stresses of the MV continue to increase until the MV is fully loaded. We note that at fully-opened state, all three cases experience similar stress level, while at just-closed and fully-loaded states, case M3 seems experiencing lower stress level compared with other two cases which may be caused by higher level compression in the leaflets (Fig. 2(h)). The transvalvular flow rates for all three cases are shown in Fig. 5. The flow rate of case M3 is higher than cases M1 and M2 in the diastolic filling phase (before 0.17 s), with M2 the lowest. When the MV begins to close, the flow rate decreases to a negative value when the closure regurgitation occurs. The regurgitation flow rates of cases M1 and M2 are similar, but much larger in case M3. Finally, the MV flow rate gradually returns to zero when the MV is fully closed. The regurgitation flow during MV closure are listed in Table 5. The results are consistent with the values from our previous work 16 . It can be found case M2 has the smallest regurgitation closing flow, while highest in case M3.
We further calculate the orifice area (OA) at fully-opened and just-closed states. In general, the higher the OA, the smaller the energy loss 55 . To calculate OA, the leaflet free-edges are first projected to the annular plane, and then the enclosed area by the project boundaries is considered to be the real OA. Table 6 shows that case M1   www.nature.com/scientificreports www.nature.com/scientificreports/ has the largest OA at fully-opened state and smallest OA at just-closed state, which indicates the impedance of case M1 in the diastolic filling is lower compared to other two cases and is closed more tightly in systole. The OA values from the three cases at fully-opened state are also within the interval (4-6 cm 2 ) given by Luo et al. 44 .
In summary, of all three models, the constitutive law M2 has the lowest leaflet displacements, lowest peak velocities at fully-opened and closed state, highest fiber stress at fully-loaded state, and lowest OA when fully-opened and the larger OA at closure. The constitutive law M1 achieves largest OA at fully-opened state  www.nature.com/scientificreports www.nature.com/scientificreports/ and lowest OA at fully-closed state, and smaller regurgitation closure flow. Therefore, the constitutive law M1 is deemed more suitable for predicting MV dynamics with FSI. In the next section, we will use the constitutive law M1 for MV leaflets to study the effects of two different material models of the chordae tendineae on MV dynamics.
The effects of different constitutive models of the chordae tendineae. Table 7 gives the displacements with different chordae constitutive laws. The average displacements are almost the same at fully-opened and fully-loaded states, while the maximum displacement from the linear model is slightly larger than that of the exponential model at fully-opened and -closed states.  Trigon 1  34  28  6  130  180  89  255  325  238   Trigon 2  18  18  12  142  213  110  280  439  260   The belly region 52  59  16  78  147  62  127  283  124   Table 3. Average stresses along fiber direction on three local regions.  Table 4. Average strains along fiber direction on three local regions.   www.nature.com/scientificreports www.nature.com/scientificreports/ We summarize the fiber stress and strain results of the exponential model in Fig. 6. Comparing with the case M1 of Figs 2 and 3, we find that the stress and strain distributions of the two chordae tendineae constitutive laws are also similar. The average stress and strain of defined three regions for the exponential constitutive law are shown in Table 8. Compared to the results for the linear model (case M1), the stress level of the exponential model is slightly higher than that of the linear model. Figure 7 shows the flow rates through the MV, although slightly higher flow rate can be achieved in the linear model during diastolic filling phase (before 0.17 s), larger regurgitation closure flow exists compared to the exponential law, suggesting that the MV closes tighter when an exponential law is chosen for chordae tendineae.  Table 7. Average and maximum displacements from two chordae tendineae models.  Table 8. Average strain and stress along fiber direction on the three local regions using an exponential model for the chordae tendineae.

Discussion
In this paper, we use the IB/FE method to study the effects of different constitutive laws on MV dynamics with fluid-structure interaction. We select three different constitutive laws of MV leaflets and two material models for chordae tendineae. Parameters of different constitutive laws of leaflets are determined by matching the biaxial stretch-stress relationship along the fiber direction and the cross-fiber direction with stretch-stress relationships derived from the constitutive law M1 16     www.nature.com/scientificreports www.nature.com/scientificreports/ regurgitation closing flow of case M1 and the exponential chordae tendineae model are lower than others. Our results may suggest that the combination of M1 constitutive law for the leaflet and the exponential law for chordae tendineae would be more suitable for simulating MV dynamics with FSI.
Through in vitro 24,25,56 and in vivo 57 studies, many researchers have also studied the mechanical properties of MV by uniaxial and biaxial stretching experiments 23,51 , and those experimental data have shown that valve leaflets and chordae tendineae have the characteristics of anisotropy and non-linearity. In papers 37,58-60 , linear elastic valve material models were used, which was impractical. Therefore, we choose fiber-reinforced constitutive laws for the MV leaflet 22,23 in this study, a common practice in current soft tissue modelling. We further find that the three selected MV leaflets constitutive laws can fit to our own porcine MV experiments very well, as shown in Fig. 10.
The strain and stress distributions along fiber direction at different time are shown in Figs 2, 3 and 6. It can be seen that for different constitutive laws, most of the leaflets regions are tensile along the fiber direction and some of regions are compressed at closed state, which is, in general, consistent with the results of the papers 22,23 . Detailed strain and stress analysis in three regions are given in Tables 3-4 and 8. In the paper 43 , the maximum fiber strain is 0.4, which is larger than our results. However, our model prediction of the circumferential stresses range seems to agree well with the measurements from in vivo measurements 25,61 . For example, the circumferential stresses 25,61 of valve leaflets range from 200 to 280 kPa at fully-loaded state. There are some differences among different constitutive laws for MV leaflets. For example, case M2 exhibits a greater stress on the posterior valves, while cases M1 and M3 experience almost the same stress level. In addition, we find that the stresses are mainly concentrated on the annulus ring and the edge of the MV leaflets during closure. When the exponential constitutive law of the chordae tendineae is used to replace the linear constitutive law, we find that the strain level of MV leaflets is larger than that of the linear model, but still within the ranges reported in 62 .
Additionally, we have compared several quantities to evaluate three different material models of MV leaflets, including the peak jet velocity, the closure regurgitation volume and the orifice area. The difference in peak velocity is minor for the three cases according to Table 2. As for the amount of regurgitation flow, case M2 has the lowest regurgitant volume and M3 has the highest regurgitant volume. Table 6 shows that case M1 has the biggest orifice area when fully opened and the smallest when just closed. The orifice area values at fully-opened state are within the interval (4-6 cm 2 ) reported 44 . 46 is employed to simulate MV dynamics in this study. The governing equations of the FSI system are

IB/FE method. The IB/FE method developed by Griffith and Luo
E E e where X = (X 1 , X 2 , X 3 ) ∈ E denotes the material (Lagrangian) coordinates in the reference configuration, x = (x 1 , x 2 , x 3 ) ∈ Ω denotes the Cartesian (Eulerian) coordinates. Ω ⊂ R 3 denotes the physical region occupied by the fluid-structure system, and E ⊂ R 3 denotes the region occupied by the immersed structure (such as the mitral valve, chordae tendineae, etc) in the reference configuration. ρ is the fluid density, p(x, t) is the Eulerian pressure, and μ is the viscosity. χ(X, t) ∈ Ω gives the physical position of material point X at time t. Therefore, the physical region occupied by the structure at time t is Ω e (t) = χ(E, t), and the physical domain occupied by the fluid at time t is Ω f (t) = Ω − Ω e (t). A three-dimensional regularised delta function δ(x) = δ(x 1 )δ(x 2 )δ(x 3 ) was used to describe the fluid-structure interaction, which implies that the IB/FE approach permits nonconforming discretization of the fluid and structure domains. = ∂Ψ ∂ P F e is the first Piola-Kirchhoff (PK) stress tensor, which is calculated from a strain-invariant based strain energy function Ψ.
The total Cauchy stress tensor of the coupled fluid-structure system is , t)) T ] is the fluid-like stress tensor. I is the identity matrix, and σ e is the elastic stress tensor, defined as is the deformation gradient and J = det(F).

Constitutive laws and parameters.
Biological tissues usually can be modeled as nonlinear elastic materials, and their material parameters could be obtained from uniaxial or biaxial tensile testing, in which tissue samples are subjected to various stretching configurations along different directions. MV anatomy experiments show that valve tissue is basically composed of fibrous tissue 51,63 , mainly collagen and elastin, and the liquid (mainly water). At low strain, the wavy structure can be extended by relatively low stress, but with the increase of strain, the fiber straightens gradually and the overall response of the structure becomes more rigid. To determine material parameters, an inverse problem is usually formulated by minimizing the differences between the predicted stretch-stress data derived from selected constitutive laws and experimentally measured data, that is where c i (i ≥ 1) are non-negative material parameters and g i are constraints of constitutive constants, i.e. >0, and σ model is calculated from some constitutive laws, σ exp are experimental measurements.
The constitutive laws of the mitral valve leaflets. In this study, three fiber-reinforced strain energy functions (cases M1, M2 and M3, Eq. 9) are chosen to characterize the mechanical responses of MV leaflets, all are based on strain invariants of I 1 and I 4 , respectively, and in which C = F T F is the Cauchy-Green deformation tensor, a 0 is the collagen fiber direction at reference state, which is an unit vector. I 1 represents the overall deformation, usually is used to describe the isotropic matrix property, and I 4 is the squared stretch along the collagen fiber direction.  where c, a, b, c c c , , 0 1 2 and c 0 , c 1 , c 2 are the non-negative material parameters. We assume collagen fibers can only bear the load when they are stretched, but not in compression, thus I 4 in Eq. 9 are replaced by = ⁎ I I max( , 1) 4 4 . These three constitutive laws of the valve leaflets are all transversely isotropic materials, case M1 was used in paper 16 , case M2 was used in paper 21 , and case M3 was from the papers 22,28 .
The corresponding Cauchy-stress tensor for the three selected strain energy functions are  in which λ is a Lagrangian multiplier to enforce the incompressibility, and B = FF T . Parameters for the strain energy function M1 are from Gao et al. 16 , which are derived from the in vitro testing on a healthy human MV carried out by Wang et al. 64 . The parameters of case M1 are shown in the Table 9. Parameters in cases M2 and M3 are determined by using a "pseudo" biaxial stretching experiments along fiber and cross-fiber directions using case M1. The fmincon function in MATLAB is employed to determine the parameters of constitutive laws M2 and M3 by using the Cauchy stress formulas (Eq. 10) by minimizing Eq. 8. The fitted stretch-stress curves for all three laws are given in , and the estimated parameters of cases M2 and M3 are listed in Table 10, including the average errors.
We further fit the three constitutive laws from Eq. 9 to an in vitro biaxial stretching experiment on porcine MV samples. In brief, fresh porcine MV samples were harvested from a local abattoir. Specimens of MV and  www.nature.com/scientificreports www.nature.com/scientificreports/ chordae were then dissected and stored in 4 °C phosphate buffer saline (PBS) before test and submerged in 37 °C PBS bath during test. Planar biaxial tensile test was conducted with a CellScale BioTester (Waterloo, ON, Canada) with 10 N load-cell on MV samples, while uniaxial tensile test (Instron Industrial Products, US) was carried out on chordae samples. The tissue specimens were stretched and released for 8 complete cycles for preconditioning until the load-displacement curve was visibly repeatable. Finally, the MV and chordae specimens were stretched to 1500 mN and 5 N to cover the physiological condition, respectively. The displacement and tensile forces were recorded and used for stress and strain analysis. The thickness of specimens was measured by a digital caliper (±0.01 mm) three times on random location before the test. We find that all three constitutive laws can fit the experimental data very well as can be seen from Fig. 10, among which cases 1, 2 and 3 represent the three sets of MV data from three hearts. The model M2 shows the best fitting although the difference is small. Corresponding R-squared values are reported in Tables 11-12, all have minor differences for the anterior and posterior leaflet respectively, which suggests that these selected three constitutive laws are suitable for characterizing MV leaflets properties. Thus, we will mainly compare their effects on the MV dynamics using a FSI solver in this study.
The constitutive laws of the chordae tendineae. Two constitutive laws are chosen for the chordae tendineae, one is the neo-hooken material model 16 , and the second one is the exponential model 28 , they are    Table 12. R-squared values of fitting the posterior leaflet with three constitutive laws in Eq. 9 to our porcine MV experiments. SSE is the residual sum of squares, SST denotes the total sum of squares. www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 11 shows the fitted stretch-stress curves with the two chordae constitutive laws to our own uniaxial tensile testing experimental data from Dr. Ma's Lab using porcine MV chordae tendineae. It can be found that nonlinear mechanical property of the chordae tendineae can only be better represented by the exponential law.
The MV model and boundary conditions. The MV model is reconstructed from a cardiac magnetic resonance (CMR) imaging of a healthy volunteer, and the leaflets are reconstructed at mid-diastole, and a pseudo-chordae structure is used because the CMR imaging cannot describe the chordal structure in vivo due to resolution limitation. Details of the MV geometry reconstruction can be found in our previous study 16,36 . Figure 12 shows the MV with the chordae tendineae, mounted in a housing and then attached to a straight tube (length: 16 cm, radius: 3.8 cm), and immersed into a fluid domain with size 10 cm × 10 cm × 16 cm, which is discretized into 80 × 80 × 128 regular grids.
An explicit version of Crank Nicolson-Adams Backward scheme is used for time stepping, which requires a relatively small time step size (10 −5 s). The IB/FE MV model is implemented within the open-source IBAMR software framework (https://github.com/IBAMR/IBAMR). The boundary conditions are the same as in paper 16 , in brief, pressure boundary conditions are applied to the inlet of the straight tube, pressure profile is shown in Fig. 13. Zero pressure boundary conditions are applied along the rest of the boundaries of the whole computational domain. The housing and the straight tube are fixed in place. CMR measured displacements of the papillary muscles are applied to the chordae attachment points where the chordae tendineae connects the papillary muscles. Further details of the MV model implementation can be found in 16 . Limitations. Finally, we mention the limitations of this study. Though we have incorporated FSI and nonlinear constitutive laws for valve leaflets and chordae tendineae. We have ignored the valve-heart interactions, which will have some impact on the dynamic loading conditions 29 . In addition, our geometric structure is based on a simplified model of chordae tendineae, whereas the realistic chordae tendineae consists of marginal, strut and basal chordae tendineae 17 . We leave investigation of these effects to future research.