Modelling and simulation of mitral valve for transapical repair applications

. Development and application of the numerical model for the simulation of human heart mitral valve (MV) transapical repair is presented. Transapical repair with neochordae implantation is a novel surgical technique allowing beating-heart correction of mitral regurgitation caused by chordae tendineae rupture through a minimally-invasive approach. In the present study, the structural ﬁnite element model of the MV decoupled from the blood ﬂow is considered. It comprises two leaﬂets and chordae tendineae described by nonlinear material model. Geometry of the model and kinematic boundary conditions for ﬁxed points of MV annulus, papillary muscles, and left ventricle apex are deﬁned by patient-speciﬁc data. Decoupled behavior of blood is speciﬁed by the time-dependent physiologic transvalvular pressure. Personalized computational modelling strategy is applied to perform virtual transapical MV repair by positioning neochordae following the real-life surgery procedure executed by surgeons. A transient analysis in time frame between end-diastole and peak systole is conducted to evaluate post-repair MV function. Computational MV simulation and modelling results provide quantitative information about the neochordae contribution to the MV function improvement and present practical value for the surgical planning of transapical MV repair.


Introduction
Mitral regurgitation (MR) is the most common heart valvular disease [23]. It is the abnormal leakage of blood backwards from the left ventricle (LV) through the mitral valve (MV) to the left atrium (LA) because of the insufficient MV closure during the cardiac cycle. This condition usually appears due to disruption of the valve restraining force structures, such as rupture of the chordae tendineae [8].
MV repair is a preferred surgical treatment for MR and has significant advantages over MV replacement, such as preservation of postoperative LV function and improved longterm survival rate [21]. However, there is no standard surgical repair technique for MR, it can depend on individual patient anatomy, pathology, and even surgeon preference [22]. Still, the use of artificial neochordae for the MV repair has become the most common method, and currently, over 40 different surgical techniques for applying neochordae to MV repair is used [11]. Expanded polytetrafluoroethylene (ePTFE) sutures are typically used to form neochordae, as their biomechanical properties are similar to those of the natural chordae [28].
Transapical repair with neochordae implantation is a novel surgical technique allowing beating-heart correction of MR through a minimally-invasive approach. NeoChord DS1000 device is passed into the beating heart through a small incision and used to implant the neochordae into the MV. The implanted neochordae are secured under proper tension to the LV apex. Echocardiographic guidance is used for both implantation and tension adjustment of the neochordae [20].
During the last three decades, a number of numerical simulation studies have shown promising outcomes while investigating the MV function [13,34,36] and evaluating novel MV surgical repair techniques [22,26,37]. Two categories of the computational models can be distinguished. Structural models use finite element method (FEM) to simulate the MV function without directly considering the blood flow, whereas fluid-structure interaction models study the interactions between the MV and the blood flow using computational fluid dynamics [18].
Numerical simulation studies of MV repair with neochordae implantation were previously published by several authors [22,26,27,32,33]. However, none of them had addressed the problem of transapical MV repair. To the best of our knowledge, our previous study [9] was the first attempt to evaluate the effects of transapical MV repair with neochordae implantation using FEM. Still, in that study, the echocardiographic data of healthy MV was used, and the insufficient MV closure was simulated by removing several chordae. Therefore, morphological comparison with the real-life surgery outcomes could not have been done.
In the present study, development of structural finite element model for the simulation of transapical MV repair is presented. Personalized computational modelling strategy is applied to perform virtual transapical repair and to evaluate post-repair MV function, comparing the simulation results with the real-life surgical procedure outcomes and with the ones previously published by other research groups.

Background
Development of the computational model for the simulation of the diseased valve function and MV repair procedure requires an understanding of its normal anatomy and function, as well as how MV is altered by the disease.

Mitral valve structure
The MV is a complex intracardiac structure that ensures unidirectional blood flow from the LA to the LV during the cardiac cycle. It consists of two leaflets (anterior (AL) and posterior (PL)) attached to the mitral annulus and anchored to the papillary muscles (PMs) by a group of thread-like branched fibrous tendons called the chordae tendineae. Both leaflets can be divided into lateral (A1, P1), middle (A2, P2), and medial (A3, P3) segments, as proposed by Carpentier [3] and widely used. However, in the present study, the modified Carpentier classification [30] is used characterized by the division of A2 and P2 segments in halves, namely middle-lateral (A2L, P2L) and middle-medial (A2M, P2M) segments (Fig. 1). This approach allows a more precise description of the location of the valve dysfunction.
The MV leaflets are attached to the mitral annulus which is located between the LA and the LV. It is a saddle-shaped fibrous ring which changes its shape throughout the cardiac cycle. The ventricular surface and the free margin of the leaflets are connected to the PMs by the branched network of the chordae tendineae. There are two PMs in the LV, each having one or multiple heads with the chordae attached. The lateral half of the leaflets is anchored by the chordae to the anterolateral PM, while the medial half has chordal attachments to the posteromedial PM. Different types of chordae can be distinguished: (a) marginal chordae are the thinnest and inserted into the free margin of the leaflets; (b) basal chordae can be found near the annular attachment of the leaflets; (c) strut chordae are particularly thick and attached to the central parts of the leaflets [17].

Mitral valve function and dysfunction
From the mechanical point of view, the MV leaflets play the role of the valve which opens and closes because of pressure differences between the LA and the LV, opening when there is higher pressure in the LA, and closing when pressure is higher in the LV. During the diastole, blood fills the LA thus increasing pressure and opening the MV. Blood begins to flow through the valve, increasing pressure in the LV and forcing the MV to close.
The chordae tendineae play an important role in the MV function, holding the leaflets in place and not allowing them to evert to the LA. This eversion is called MV prolapse. Specifically, basal and strut chordae provide only basic structural support and are more involved in the LV function, whereas marginal chordae prevent prolapse of the leaflet margin [24]. For this reason, ruptured basal and strut chordae might not compromise the coaptation of the leaflets. However, rupture of marginal chordae almost always results in loss of normal coaptation and valvular leakage [5].

Materials and methods
The modelling strategy applied for virtual transapical MV repair consisted of following steps: (i) 3D image data acquisition before and after real-life surgery; (ii) manual tracing of the MV structure on acquired images and reconstruction of the 3D MV geometry; (iii) mesh generation and implementation of mechanical properties, initial and boundary conditions; (iv) virtual transapical MV repair following the real-life surgery procedure; (v) transient simulation in Abaqus/Explicit.

Surgery and data acquisition
The data for the present study was obtained during transapical MV repair surgery performed at Vilnius University Hospital Santaros Klinikos (Lithuania) using NeoChord DS1000 device. This device is designed to deliver neochordae under beating-heart conditions with transesophageal echocardiography (TEE) guidance. During TEE, the ultrasonic transducer is placed into the esophagus of the subject to get various real-time views of the cardiac structures.
The patient was diagnosed with MV prolapse of the P2L segment and moderate to severe MR (grade 3) due to the rupture of marginal chordae of the posterior leaflet. During the surgery, a total of three neochordae were implanted into the posterior leaflet and secured to the LV apex: two neochordae were inserted into the P2L segment and one into the P2M segment for better support. ePTFE CV-4 type sutures with a diameter of 0.307 mm were used to form the neochordae. Unfortunately, criss-crossing between neochordae and native chordae of the anterior leaflet was noticed, and in a few minutes, rupture of native chordae occurred causing prolapse of the A2L segment. For this reason, three more neochordae were implanted into the anterior leaflet. The implantation of neochordae has ensured decent support for both posterior and anterior leaflets, as well as reduced MR to mild (grade 1).
During the surgery, data sets of the movement of pre-and post-operative MV were obtained using Philips EPIQ 7C ultrasound machine and stored as a volumetric medical image (VolDICOM) for further processing.

Mitral valve geometry
The reconstruction of the prolapsing MV geometry was performed using the custom platform developed in MATLAB by Biomechanics Research Group of Politecnico di Milano (Italy) [31]. The Cartesian reference system was positioned within the acquired volumetric image data such that z-axis would go through the LV apex and the MV centroid, while xand y-axes would be perpendicular to the MV annular plane. The 3D data was cropped into 18 radial planes evenly distributed around z-axis. For each leaflet, the positions of the annulus and the free margin were manually identified on every radial plane. In addition, multiple points were selected to describe the profile of the MV leaflets (Fig. 2). The identified points were interpolated with a cubic spline of 32 points, uniformly distributed along the length of each leaflet ( Fig. 3(a)). This procedure was performed on the enddiastole frame of TEE image data, chosen as the last frame before the MV starts to close, with a hypothesis that at this point in time the MV is approximately unloaded.
In order to avoid artifacts due to manual tracing, Fourier approximation was used. A cylindrical reference system was set with the origin in the centroid of the manually identified annulus points. In this local reference system, the radial (ρ) and axial (z) positions of every spline point were approximated through 4th order Fourier functions of the angular position (φ): where a 0 , a k , b k , c 0 , c k , d k are the coefficients of the resulting Fourier function. Every Fourier function was upsampled to 200 points, and a 3D point cloud was generated, consisting of 32 levels running from the annulus to the free margin, each with 200 points, uniformly distributed along the circumference of the MV (Fig. 3(b)). Finally, the coordinates of all identified points were transformed back into Cartesian coordinates with the origin and z-axis coinciding with those of the cylindrical reference system. http://www.journals.vu.lt/nonlinear-analysis

Mechanical properties of the leaflets
The MV leaflets were modelled as a thin membrane described by nonlinear shell elements. The plane stress assumption was made meaning that Cauchy stress tensor σ component σ 33 is equal to zero.
Generated point cloud was connected into a mesh of 3-node finite-strain shell elements (type S3R in Abaqus). 18 900 shell elements were created. Regionally varying thickness was assigned to the leaflets, as suggested by Kunzelman et al. [14], with an average value of 1.32 mm for the anterior and 1.26 mm for the posterior leaflet. Density for the MV leaflets was assumed to be 10.4 g/cm 3 , which is ten times higher than the real density, thus taking into account the inertial effects of the blood volume spanned and moved by the MV during its closure. Such volume is about ten times larger than the volume of the MV leaflets.
Structurally, the MV leaflets present a soft tissue abundant in collagen fibers oriented parallel to the annulus thus increasing stiffness in this direction. Therefore, their mechanical properties were assumed as nonlinear and anisotropic. Since there are several proposed constitutive models which represent such behavior of the MV leaflets [19,25,38], the one suggested by Lee et al. [16] was chosen. It considers the effects of material anisotropy at the tissue level enabling the investigation of how the material anisotropy at such level affects the mitral valve functional response at the organ level. Therefore, the stress-strain behavior of the leaflets was described using hyperelastic strain energy density function: where I 1 = tr(C) and I 4 = NCN, in which C = F T F is the right Cauchy-Green strain tensor, N is the unit vector identifying collagen fibers direction in the unloaded configuration, F = ∂x/∂X is the deformation gradient tensor defined as the derivative of the position in the deformed configuration x with respect to the position in the original (undeformed) configuration X. β is a constitutive material parameter ranging from 0 to 1 and related to the level of material anisotropy, and c 0 , c 1 , c 2 , c 3 are the remaining material parameters. The constitutive material parameters (Table 1) were obtained by stress-strain curve fitting to biaxial test data reported by May-Newman and Yin [19] using a nonlinear least square method.
The constitutive relation (3) was implemented in a user-defined VUMAT subroutine within Abaqus in which the derivatives of strain energy Ψ with respect to the invariants I 1 and I 4 are evaluated: and the components of the Cauchy stress tensor σ are returned as an output: where C 11 , C 12 , C 22 are the components of the right Cauchy-Green strain tensor C, λ 11 is the stretch ratio in the direction parallel to the MV annulus, λ 12 is the shear stretch ratio, and λ 22 is the stretch ratio in the direction perpendicular to the MV annulus.

Mechanical properties of the chordae tendineae
The MV model was completed by adding the PMs and the chordae tendineae (Fig.3(c)). The PMs were assumed as single nodes neglecting the presence of multiple heads in the same PM. Therefore, the coordinates of only PM tips were manually identified on the end-diastole frame of TEE image data. A branched network of marginal, basal, and strut chordae was created connecting the PM tips and the leaflets. The insertion points on the leaflets were determined in accordance with ex vivo findings [15] and indications by clinicians [31]. The chordae tendineae were modelled as 2-node tension-only truss elements (type T3D2 in Abaqus) and assigned constant cross-sectional area values of 0.40, 0.79, and 1.15 mm 2 for marginal, basal, and strut chordae, respectively [35]. Density for the chordae tendineae was assumed the same as for the MV leaflets and equal to 10.4 g/cm 3 . 247 truss elements were created. The chordae tendineae are mainly composed of collagen fibers aligned along the longitude of these tendons thus making them both stiff and robust. Therefore, mechanical properties of the chordae were assumed as nonlinear and isotropic. Since the chordae tendineae undergo only small or moderately large strains (< 100%) during the cardiac cycle, their mechanical behavior was described using 2nd order polynomial hyperelastic model from Abaqus material model library. The constitutive material parameters were automatically computed by Abaqus from uniaxial test data reported by Kunzelman and Cochran [12].
PM tips were treated as single nodes without physical properties.  of TEE image data before the MV starts to close. At this point, in time, the MV can be assumed as approximately unloaded. Therefore, no prestress was included in the model. The following boundary conditions were applied to the MV model:

Initial and boundary conditions
(i) Patient-specific kinematic boundary conditions were applied to achieve the personalization of the presented modelling strategy. The motion of the MV annulus and the PMs were traced up to the peak systole frame of TEE image data. At peak systole, ventricular volume decreases to the lowest level, and ventricular pressure reaches its maximum, while the MV remains closed. The annular contraction and the movement of the PMs were modelled by applying nodal displacements to the points of the annulus and the PM tips ( Fig. 4(a)). (ii) Decoupled behavior of blood was specified by applying the time-dependent physiologic transvalvular pressure ranging between 0 and 16 kPa on the ventricular surface of the leaflets (Fig. 4(b)). (iii) To capture the contact of the leaflets, a general contact algorithm available in Abaqus with a scale penalty method and a friction coefficient of 0.05 was set between the atrial sides of the leaflets.

Virtual repair
The chordae tendineae were removed from the edge of the P2L segment of the created model to replicate MV prolapse seen on pre-operative TEE images. During transapical MV repair, the neochordae implanted into the leaflets are secured to the LV apex. For this reason, the position of the LV apex as a single node was manually identified on the end-diastole frame of TEE image data. Following the real-life surgery, three virtual neochordae were implanted into the posterior leaflet, connecting the LV apex and the free margin of the leaflet. The bulging height of the posterior leaflet from the annular plane at peak systole was identified. The lengths of the neochordae were determined by subtracting this height from a distance between the LV apex and the free margin at peak systole (Fig. 5). The motion of the LV apex was traced up to the peak systole frame of TEE image data and applied to the LV apex as a nodal displacement. After the implantation of the neochordae into the posterior leaflet, the same procedure was repeated with the A2L segment thus allowing to recreate prolapse of the anterior leaflet which occurred during the surgery due to the rupture of the native chordae. The neochordae were modelled as truss elements (type T3D2 in Abaqus) and assigned constant cross-sectional area value of 0.074 mm 2 , which corresponds to the crosssectional area of ePTFE CV-4 type suture, typically used for transapical MV repair. Density for the neochordae was assumed to be ten times higher than the real density and equal 22.0 g/cm 3 , thus taking into account the inertial effects of the blood. Mechanical properties of the sutures were described as nonlinear and isotropic using 2nd order polynomial hyperelastic model. The constitutive material parameters c ij were computed from uniaxial test data reported by Dang et al. [6].
The LV apex was modelled as a single node without physical properties.

Simulation setup
Four finite element models were prepared for the transient simulation in Abaqus representing: (i) prolapse of the posterior leaflet (PLP); (ii) virtual repair with neochordae implantation into the posterior leaflet (PLVR); (ii) prolapse of the anterior leaflet which occurred during the surgery (ALP); and (iv) virtual repair with neochordae implantation into the anterior leaflet (ALVR). Abaqus/Explicit was chosen to perform the simulations since it employs central-difference integration scheme to solve nonlinear systems with complex contacts without iterating and requires much less memory and disk space compared to implicit integration, therefore, is more suitable for solving complex dynamic problems for nonlinear materials. For each model, the input file for the analysis in Abaqus/Explicit was prepared. The largest time increment ∆t for the integration procedure to remain stable was automatically calculated by global estimation algorithm in Abaqus to be equal 2.5·10 −6 s. However, this estimated ∆t is only approximate and in most cases actual stable time increment should be less than this. For this reason, automatic time incrementation with the maximum allowed time increment of 0.5 · 10 −6 s was specified for all simulations.
Simulations were carried out one by one thus allowing to recreate the real-life surgery process.

Morphological compliance
In the present study, a rare bileaflet transapical MV repair which occurs in less than 5% of cases [4] was investigated. Personalized computational modelling strategy was applied to perform virtual transapical MV repair. The simulation results were used to evaluate morphological compliance between the virtual models and the TEE image data. This was done by visually comparing the computed systolic leaflet configurations of virtual PLP and ALVR models with the pre-and post-operative echocardiographic images. Figure 6 shows a good agreement in deformed leaflet morphology between the simulation results and the echocardiographic data.

Computational results
After each simulation (i.e. prolapsing posterior leaflet followed by virtual repair and prolapsing anterior leaflet followed by virtual repair), the following parameters were identified at peak systole to evaluate the effects of virtual transapical MV repair: (i) maximum displacement along z-axis of the MV free margin (u z,max ); (ii) coaptation area (CoA), defined as the contact area of the leaflets, and coaptation length (CoL), defined as the length of the leaflet apposition in the prolapsing region; (iii) tension forces in native chordae (F CT ) and neochordae (F NC ); and (iv) peak value of the maximum principal stress acting on the leaflets (S 1,peak ). The values of the computed parameters are shown in Table 2.  Both PLVR and ALVR repositioned prolapsing segments of the leaflets below the annular plane, resulting in the elimination of MV prolapse. In PLP model, the maximum z-axis displacement of 18.6 mm appeared on the free margin of the prolapsing P2L segment. Following PLVR, prolapse of the posterior leaflet was eliminated, and the maximum displacement value was reduced to 13.5 mm. In ALP model, prolapse of the anterior leaflet occurred, and the maximum z-axis displacement of 19.1 mm was noticed on the free margin of the A2L segment. ALVR eliminated prolapse and reduced the maximum displacement value to 12.9 mm.

Coaptation area and coaptation length
In the initial model, a lack of contact between A2L and P2L segments due to the prolapsing posterior leaflet was present. After PLVR, a contact between the segments was recovered restoring the coaptation length of 5.8 mm and increasing the coaptation area by 21.5%. However, after the rupture of the A2L native chordae, the bulging of the anterior leaflet occurred resulting in a decrease of the coaptation area by 15.9%. A good contact between the leaflets was restored after ALVR increasing the coaptation length from none to 6.4 mm. Compared to the PPL model, a total increase of the coaptation area by 41.9% was achieved during both virtual repair procedures.

Tension forces
The maximum tension forces in the chordae tendineae were found in the marginal chordae neighboring the prolapsing segment. In PLP model, the peak value of tension force in the marginal chorda inserted into the border of P1-P2L segments was 1.1 N, whereas in ALP model the maximum tension force of 1.08 N occurred in the marginal chorda inserted into the edge of A1-A2L segments. After each virtual repair procedure, these values were reduced to 0.65 and 0.5 N, respectively. Tension forces in the neochordae were calculated after PLVR and ALVR. Following PLVR, the maximum value of the force exerted by the neochordae implanted into the posterior leaflet was 0.96 N. This value was decreased to 0.86 N after ALVR, whereas the maximum tension force in the neochordae implanted into the anterior leaflet was calculated to be 0.88 N.
The change of tension in the native chordae and neochordae as a whole can be represented as an alteration of reaction forces acting on the PMs and LV apex during transient analysis. Figure 7 shows that virtual repair procedures reduced reaction forces on the PMs since tension in the native chordae was partially transferred to the neochordae, resulting in the appearance of reaction force on the LV apex.

Stress analysis
Stress distribution across the leaflets of the prolapsing and virtually repaired MV models are shown in Figs. 8, 9. All simulations specified the peak value of the maximum principal stress to be on the MV free margin. In PLP model, the high stress concentration area occurred near the prolapsing P2L segment with a maximum value of 1.09 MPa. Following virtual repair, this high stress concentration was significantly reduced. However, the increase of the stress near the attachment points of the neochordae was noticed with a peak value of 0.71 MPa. After the rupture of the A2L chordae, new high stress concentration areas appeared on the anterior leaflet. As previously, the highest stress values occurred near the prolapsing http://www.journals.vu.lt/nonlinear-analysis ALP ALVR Figure 9. Maximum principal stress distribution across the mitral valve leaflets before and after implantation of the neochordae into the anterior leaflet. Areas with stress peak values are circled in red.
segment along the free margin of the leaflet with a maximum value of 0.85 MPa. In addition, stress increase was noticed in the middle of the bulging region of the leaflet. Following ALVR, these stress values were reduced, still, as before, the increase of the stress appeared near the neochordal attachment points with a peak value of 0.79 MPa.

Discussion
The parameters calculated during transient simulations were compared to those found in the literature. Numerical simulation studies quantitatively investigating different neochordal implantation techniques were used for the comparison. The maximum displacement along z-axis of the MV free margin and the coaptation area of the leaflets are the parameters strongly dependent on the geometry of the MV. Therefore, these parameters are hardly comparable with the ones found in the literature. Still, the post-repair maximum z-axis displacement and coaptation area recovery values are similar to the ones reported by Sturla et al. [32]. In addition, both virtual repair procedures restored a sufficient coaptation length larger than 5 mm which is considered to be a minimum coaptation length required to ensure adequate MV function [7].
In terms of tension forces in the native chordae, the calculated values are slightly higher than the ones reported by Rim et al. [27] though they do not exceed the chordal failure tension values found in the literature [29]. This difference could be due to the higher systolic pressure used in the present study with a peak value of 16 kPa, compared to 12 kPa used by Rim and colleagues. Regarding the tension forces in the neochordae, the calculated values are comparable with the ones measured in vitro [2] and in vivo [1] by Bajona et al. Tension forces, estimated in the current study, fall into the reported range between 0.6 and 1.1 N.
The computed stress patterns are similar to the ones reported by the other authors [27,32]. In both models with MV prolapse, high stress concentration areas occurred near the prolapsing segments, whereas following virtual repair procedures the increase of the stress near the attachment points of the neochordae appeared. The calculated peak values of the maximum principal stress did not exceed the previously reported failure stress values [10].
The comparison to the existing results shows that all the computed parameters have a good agreement with those reported by the other authors and do not exceed any critical values found in the literature.

Conclusions
In the present study, personalized computational modelling strategy for transapical MV repair application allowed to determine quantitative information about the contribution of the neochordae to the MV function improvement. Virtually repaired models demonstrated eliminated prolapse, improved leaflet coaptation, decreased chordal tension forces, and reduced excessive leaflet stress while maintaining the morphological compliance with the echocardiographic data.
The obtained computational results indicate the possibility of using the presented modelling strategy as a tool for predicting post-operative MV morphology and coaptation of the leaflets as the main indicators of MV prolapse elimination. Therefore, it shows the potential to have practical value in scenarios of clinical relevance.