|Home | About | Journals | Submit | Contact Us | Français|
Numerical models of native heart valves are being used to study valve biomechanics to aid design and development of repair procedures and replacement devices. These models have evolved from simple two-dimensional approximations to complex three-dimensional, fully coupled fluid-structure interaction (FSI) systems. Such simulations are useful for predicting the mechanical and hemodynamic loading on implanted valve devices. A current challenge for improving the accuracy of these predictions is choosing and implementing modeling boundary conditions. In order to address this challenge, we are utilizing an advanced in-vitro system to validate FSI conditions for the mitral valve system. Explanted ovine mitral valves were mounted in an in vitro setup, and structural data for the mitral valve was acquired with μCT. Experimental data from the in-vitro ovine mitral valve system were used to validate the computational model. As the valve closes, the hemodynamic data, high speed leaflet dynamics, and force vectors from the in-vitro system were compared to the results of the FSI simulation computational model. The total force of 2.6 N per papillary muscle is matched by the computational model. In vitro and in vivo force measurements enable validating and adjusting material parameters to improve the accuracy of computational models. The simulations can then be used to answer questions that are otherwise not possible to investigate experimentally. This work is important to maximize the validity of computational models of not just the mitral valve, but any biomechanical aspect using computational simulation in designing medical devices.
Beginning with Da Vinci, the structure and function of the mitral valve have been of great interest to anatomists, physicians, surgeons, scientists, and engineers. Yet, a complete understanding of the complex mechanical function of the normal mitral valve remains elusive. The anatomy of the mitral valve (MV) includes two asymmetric leaflets attached to the left ventricle at the annulus, and numerous chordae tendineae that serve to attach the leaflets to the left ventricle via the papillary muscles. The chordae connect to both the free edge of the leaflet (primary chords) but also to the ventricular surface (secondary chords) (Fig. 1). Some smaller chords (tertiary) do not connect to the leaflet at all but rather provide connections between adjacent chords. The function of the valve throughout the cardiac cycle is complex, and involves interaction between all of the aforementioned structures.
Common diseases of the MV are related to dysfunction the valvular apparatus as a whole, or due to alteration or disruption of particular parts of the valve. Stenosis (resistance to forward blood flow), regurgitation (reversal of flow), and prolapse (displacement of leaflets towards the atrium) are pathologic conditions that can result from valvular tissue abnormalities and disruption. To treat such pathology, MV repair is the most popular and most reliable surgical treatment. However, one of the unsolved problems in MV repair is the ability to predict the optimal repair strategy for each patient. Experimental studies have helped to improve repair techniques, and now, computational simulations are becoming increasingly utilized in an attempt to understand the complex MV dynamics.
In the case of our group’s ongoing analyses of the mitral valve, finite element analysis was introduced three decades ago to begin to understand the mechanical stresses in the valve in normal, diseased,[2,3,4] and surgically repaired states.[5,6,7,8] More recently, we have utilized fluid-structure interaction (FSI) analysis to understand the complex dynamics mitral valvular function.[9,10,11,12] Of course, many other research groups have since also utilized finite element analysis, computational fluid dynamics, and fluid-structure interaction model in order to understand various aspects of mitral valve function (see for a comprehensive review).
Most of these models utilized MV geometries obtained either by parametric modeling, by image reconstruction of markers placed on the valvular apparatus, or from excised valve specimens. More recently, MV modeling strategies have been transitioning to morphologically realistic patient-specific MV modeling utilizing real-time 3D imaging data.[14,15,16] However, one of the obstacles preventing accurate patient specific modeling, has been the inability to accurately represent the complex chordal structure of the native valve.
For the sake of efficiency, current computational models of the MV often use simplified geometry. As current imaging modalities cannot provide a complete description of structure and distribution of the chordae in vivo, general anatomic information from clinical studies is commonly implemented into computational MV modeling for the chordae. The models commonly take advantage of the chordal beam-type structure to justify the use of beam elements (i.e. have no bending degrees of freedom) and the thinness of the leaflets to justify the use of shell elements.[17,18,12,19,20,21,22,23,24,25] Some models do represent the leaflets with 3D elements (e.g. tetrahedral) but often only utilizing a single layer in order to reduce computational time. Further, the detailed 3-D branching nature of the chords has not been well-represented. All such assumptions sacrifice the accuracy of the computational results. The chordal structure is an important functional part of the mitral valve and should not be omitted or simplified, especially when investigation of the pathologies related to chordae tendineae is of interest. Improvement in computer technology has been exponential in recent decades, and now allows the use of more detailed models, without unnecessary limitations in regard to chordal geometry and its complexity.
The purpose of this investigation was to develop an anatomically realistic FSI model of the mitral valve that includes leaflet and chordal elements that accurately recreate the natural anatomy of the mitral valve, based on CT images of a normal valve. In addition, we set out to include the annular and papillary muscle components in order to accurately simulate physiologic conditions. We examined the results with a particular focus on determining the resulting forces on the chordae tendineae throughout the cardiac cycle. Lastly, the computational studies were compared with experimental in-vitro measurements, to validate our FSI model results. In order to fully validate a computational model, as the valve closes, the hemodynamic data, leaflet dynamics and force vectors can be compared between the in-vitro experiments and the FSI model. This paper focuses on the papillary muscle forces comparing their magnitude and direction.
It is anticipated that further advancements of clinical imaging modalities, coupled with the next generation of computational techniques will enable physiologically more realistic simulations.
To develop the “valve-specific” geometry for the FSI simulation of the mitral valve with preserved chordal structure, we utilized an ovine mitral valve that was obtained from Superior Farms California. The process of model development included 1) CT scanning, 2) image processing to develop 3D model, and 3) high quality, robust mesh generation (Fig. 2). We then utilized an FSI approach to evaluate opening and closing of the valve, and to determine the time dependent chordal forces. To validate the results, explanted ovine mitral valves were mounted in an in vitro setup and subjected to normal cardiac cycle pressures.
The geometric data for the MV model was acquired with micro-computed tomography (μCT). The atrial chamber and aortic section were removed from the left heart simulator, and the left ventricle was fixtured to the μCT gantry using a custom adaptor plate. The mitral valve geometry was retained from the earlier hemodynamic characterization. The mitral valve was scanned in air under ambient pressure (~1 atm, mitral leaflets were open) using a vivaCT 40 system (Scanco Medical AG; Brüttisellen, Switzerland). The geometry was acquired at 39 μm voxel size (~600 image slices) using scanning parameters optimized for low density soft tissues (55 keV energy, 109 μA intensity, and 300 ms integration time).
A simple threshold was applied to acquiredμCT images to segment the image into in side, outside and disputed regions. Given those region assignments, a PowerWatershed algorithm[27,28] was applied to perform the final segmentation. Subsequently, Marching Cubes was applied to the binary image to extract a triangulated surface. All image processing and surface reconstruction was performed in the MATLAB/Octave-based interoperative toolkit BioGeom.
Following surface extraction, ten iterations of volume-conserving smoothing  was applied to the surface to remove staircasing while conserving geometry and volume. Subsequently, the surface was remeshed such that surface triangles were proportional to the gradient-limited feature size (GLFS)[31,32,33] The resulting geometry was segmented into anterior and poster leaflet section, anterolateral and posteromedial papillary muscles, and individual branched chordae tendineae (Fig. 3). Segmentation of the surface mesh allows for the prescription of anatomically specific material parameters. Finally, a layered, scale-invariant tetrahedral grid was generated such that there were three layers of tetrahedra through the thickness of the anterior and posterior leaflets.[31,32,33] The final mesh consisted of 217,462 tetrahedra. At run time, the the leaflet tetrahedra were converted to quadratic elements by the IMPETUS solver.
The mitral leaflet tissue material behavior (Equation 2) is characterized as an oriented entangled population of crimped collagen fibers embedded in an isotropic phase of gel-like glycoproteins and a network of isotropic elastin (Fig. 4).[11,10]
Determination of collagen fiber orientation in the mitral valve is both critical and challenging. Currently, no perfect method exists for the comprehensive experimental measurement of the three dimensional collagen architecture that includes all of the sub-structures of the mitral valve. Two dimensional characterization that focuses on the leaflets unfolded into a plane has been carried out with small angle light scattering (SALS).[34,35] Recent efforts have been made to map 2D SALS data to 3D image-derived models, using control points and interpolation. However, these mappings exclude the mitral valve chordae and artefacts from excision such as the relaxation of collagen fibers and subsequent re-orientation remain. An added challenge is how to determine fiber orientations non-invasively.
The approach we have adopted is best described as a geometric technique that is subject to experimentally determined boundary conditions. It too is not perfect but it does have some advantages: it is fast, it is fundamentally three dimensional, it requires no invasive experimental techniques and it accommodates a realistic transition between chordae and leaflets.
The gist of the method is that we have some certainty with respect to some areas of alignment. These are our boundary conditions. For example, we know that the fiber orientation is mostly aligned with the axis of the chordae, and the free edges of the leaflets. We also know that fiber directions will be approximately in the plane of the leaflets. From our previous SALS studies of fiber orientation, we know that fibers will be oriented approximately vertically at the fibrous trigones, and approximately circumferentially along the annulus at the midpoints of the anterior and posterior leaflets. Lastly, we assume that to be mechanically optimal the fiber field must be slowly varying between these fixed directions. Thus, we solve a modified Laplace problem to determine the unspecified fiber directions in between.
More specifically, let Ω be the volume over which fiber orientations need be determined, and S be the set of seed points where the fiber directions are known. We solve Laplace’s equation for the 3×3 matrix M over Ω with Dirichlet boundary data constructed from :
Note that we must solve for rather than the fiber direction itself to negate the confounding effect of polarity. In other words, we wish, for example, to treat < 1, 0, 0 > as the same direction as < −1, 0, 0 >. The solution M(x) over Ω \ S is then decomposed into an eigenvalue decomposition:
with λ1 ≥ λ2 ≥ λ3 ≥ 0. The first (orthonormal) column vector e1 in E is taken to be the fiber direction at each point x in Ω \ S, where the eigenvalue solution simply unpacks the square matrix M.
The constitutive model for the material properties utilized a three-dimensional splay invariant, based on an approximation of a three-dimensional Gaussian distribution of fibers. The model is defined in terms of the 2nd Piola-Kirchhoff stress, S, as follows
where J is Jacobian, κ is bulk modulus, σi is passive fiber stress-strain rule for the ith (1 or 2) fiber population, ϵ is active fiber stress-strain rule, Dev is deviatoric projection operator, and C is right-Cauchy deformation, and λi is dispersed fourth invariant defined as
where is the isochoric part of the right-Cauchy deformation, K is called the dispersion tensor or anisotropy tensor and is given in global coordinates; note, as K embodies the concept of dispersion, λ is not a stretch in a classical sense. The matrix K is governed by splay parameters described elsewhere. The biaxial stress-strain behavior of the mitral leaflets (Fig. 6) is characterized by the structural constitutive equation.
The passive fiber model is defined in the fiber coordinate system, and is given by a micro-structural model based on the physiology of crimped collagen fibers (see Alg. 1); specifically, given a fiber stretch λ, this model returns the true stress σ/λ and tangent modulus dσ/dλ of the fiber. There are three physiologic parameters (material constants defined at the top of the algorithm) that the user must supply. The discontinuity in the dσ/dλ curve indicates that the model presented in Alg. 1 predicts a stress σ response that is continuous and once-differentiable in stretch λ up to fiber failure.
Specific parameter values used in this model are given below. For all structures, the bulk modulus, κ, was 2.804 kPa; the isotropic parameter μ was 10.0 kPa and the fiber stiffness, Ef was 4.84 kPa. Dimensionless parameters H0/r0 and R0/r0 were 14.5 and 1.75 for the anterior leaflet, 7.0 and 1.0 for the posterior leaflet, and 8900.0 and 25.0 for the papillary muscles. For all chordae, R0/r0 was 5.0. H0/r0 varied for most chordae, ranging from 1.35 for the strut chordae to 90 for the marginal chords.
A major challenge in modeling the fluid-structure interaction of cardiac valves is the dealing with the contact between the leaflets at coaptation. Various approaches have been adopted to address this challenge, including interface capturing approaches and boundary tracking approaches. In this study, we adopted an approach based on smooth particle hydrodynamics, wherein the fluid is modeled as a collection of discrete particles.
Specifically, at run time 346,676 discrete particles were created in the fluid domain. Fluid motion and boundary interaction was solved with the the IMPETUS Afea SPH Solver®, while large deformation in the solid mitral valve was simultaneously solved with the IMPETUS Afea Solver®. All solid elements were fully integrated removing the possibility of hourglass modes and element inversion that plagues the classic under-integrated elements.
Both fluid and solid domains and their interaction were solved with an explicit integration scheme. All simulations were solved on a standard workstation. Parallel acceleration was achieved with a Tesla K40 GPU with 12 GB of Graphic DDR memory and 2880 CUDA Cores. Both FE and SPH solvers use the GPU for parallel processing thus avoiding the need for a large compute cluster. To confirm that convergence was reached, h-refinement of the finite element mesh was performed and the solution was found to yield same results.
The fluid particles were confined in a pipe-like rigid structure surrounding the MV model and prescribed velocity boundary conditions were applied to the open ends to develop clinically relevant pressures, see Fig. 8. The prescribed velocity profiles used in the simulations are based on the data used in the in-vitro measurements and shown in Fig. 11.
The points at the bottom of the papillary muscles were fixed in all three directions, as well as the points on the annular attachment.
The leaflet and chordal stress computed by the FSI simulations showed a distinct increase in stress magnitude at full closure (Fig. 7b), as compared to the open state (Fig. 7a). The stress was not symmetrically distributed, a result of the valve specific geometry, which is itself not symmetric. The local maxima of the principal stress values (1 MPa) were found primarily on the chordae tendineae and their connections to the leaflets.
Fig. 8 shows the displacement of the fluid particles in z direction at the peak systole where the closure of the MV is reached, T=Tsys. It can be seen that, at this time point, the pistons driving the flow of the particles have moved only by 6 mm for the MV to reach the closure. Thus, not many particles flowed backward to the left atrium before the MV closed. Also, the least moving particles can be found on the outside of the valve leaflets and chordae, where the pressure develops and causes the valve to close.
To evaluate the chordal force results, we examined the forces at the individual attachment points of the chords to each papillary muscle. The maximum force was found to be 0.80 N at the attachment of a strut chord attached to the anterior papillary muscle (Fig. 9).
Further, we examined the forces at each chordal attachment point at the beginning and end of systole for both the APM and PPM (Fig. 10). It can be observed that despite significant changes in the directions of the individual vectors (red arrows), the overall averaged directions (blue arrows), did not vary significantly.
To validate our results, we compared the computational FSI results to the experimental results from the in-vitro system (Fig. 11) at a transvalvular pressure of 100 mmHg. There was excellent agreement both in terms of direction and magnitude of the overall force vectors for the anterior and posterior papillary muscle attachments (Fig. 12). The total size scaled force magnitude was 2.6 N per papillary muscle for both the FSI and in-vitro results, as previously also found in other studies.[11,25]
Further, the time course chordal attachment forces along all three axes were in very good agreement between computational FSI and experimental result (Fig. 13). The total force closely follows the left ventricular pressure curve (Fig. 14) with a maximum of 2.6 N at a transvalvular pressure of 100 mmHg, as noted previously. It is also worth noticing the the pressure values start rising before the force values, which is physiologically correct.
Ventricular dilatation causes a displacement of the PM tips away from the annulus in the apical, lateral, and posterior directions, thereby inflicting an increased tethering force on the chordae tendineae.[39,40,41,42]
In this paper, we compared two different valves from two different species to validate the computational results. We used an ovine mitral valve to acquire the structural data for the model and compared it with PM force measurements performed in vitro on a porcine mitral valve. The gross structures of mitral valves across species, especially human, porcine, ovine, are similar, which is why the porcine and ovine hearts are excellent models for human cardiac anatomy. Independent of species, the mitral valves have two main leaflets, commissures, and a network of primary, secondary and tertiary chordae tendineae anchoring to the papillary muscles which are positioned at about the same species-independent location. We have shown that when the valve closes, the resulting (averaged) forces point in the same direction and, scaled to size, are about the same magnitude. Therefore, we may conclude that although the porcine, ovine and human valves are different in various subtle ways, in the large scale comparison between the biomechanics of porcine and ovine valve (and hence human), it does not have significant impact on the averaged PM force directions and magnitudes. This approach of analyzing papillary muscle mechanics for a ‘PORVINE’ heart, i.e. analysis created by superimposing PORcine and oVINE data, has been previously suggested by Ingels and Karlsson.
In vitro and in vivo force measurements are important in validating and adjusting material parameters in computational models. The models will answer questions that are otherwise not possible to investigate experimentally, for example, when the force is measured to decrease in one part of the valve, it must be increased in another part to meet static equilibrium. This work is important to maximize the validity of computational models of the mitral valve. Extracting patient specific geometries of MV while preserving its chordal structure is challenging and recent computational models therefore simplify the models used. Some experimental research work, such as by Jensen et al., can only be used for validation of computational models and algorithms if a subject-specific comprehensive MV model with 3D chordal structure is utilized. Computational models with simplified chordal structures would not provide the information necessary for comparison with the experimental data.
To mention the limitations of this methods paper, the method presented is used with and results are shown based on single simulation/experiment, i.e. this is a case study. Thus, although accurate outcome is shown, for statistically meaningful results more simulations are necessary. Moreover, even though the model used is very detailed as is, higher level of detail could still be reached. However, for this purpose, where the forces are measured at the papillary muscle tips, it is not expected to have significant impact. If some detail is lost in the image processing phase, it is largely limited to the area of individual attachment points of the chords to leaflets, not papillary muscles. Furthermore, even though the parameters for the material model are provided above, considering the applied pressure on the leaflets and the forces in chordae counteracting this to have equilibrium, this balance can be achieved for a large range of material parameters employed for the leaflets and chordae. Measuring force and displacement on separate tissue parts such as individual chordae and leaflets, however, enables improved material parameters and modeling.
This study was supported by a grant from the National Heart Lung and Blood Institute (R01-HL092926).
Given H0/r0, R0/r0, Ef, and λu, where H0 is the initial wavelength of crimp, R0 is the initial amplitude of crimp, r0 is the initial fibril radius, Ef is the elastic modulus of the fiber in the linear region, and λu is its ultimate stretch, then:
Compute the constant parameters:
where L0 is the chord length of helix over one wavelength, while Λ is the stretch and Es is the secant modulus at the transition between the toe and linear regions. If λ ≤ Λ Then
Else If Λ ≤ λ ≤ λu Then
Else Fibril Failure
Return σ/λ and dσ/dλ.
Conflict of interest
No benefits in any form have been or will be received from a commercial party related directly or indirectly to the subject of this manuscript.