Fast Image-Based Model of Mitral Valve Closure for Surgical Planning

,


Introduction
The mitral valve is the most complex of the four heart valves and is the one most often associated with disease [21]. It consists of two leaflets that open and close as the heart beats to ensure one-way flow of blood into the left ventricle. The leaflets are restrained by fibrous chords during closure. See Figure 1. Mitral regurgitation (MR) occurs when the valve fails to close adequately during ventricular contraction and blood leaks backward through the incompetent valve. It can be caused by ischemic heart disease, dilated cardiomyopathy, rheumatic valve disease, or infection [5]. MR can lead to heart failure if left untreated, and the only effective treatment is surgery. The two primary surgical treatment options are repair of the native mitral valve tissue and replacement with a prosthetic valve. Repair has been shown to result in better function and long-term survival than replacement [9,15,20], so surgical repair of the mitral valve is preferable to valve replacement for the majority of patients who require treatment for MR [7]. However, replacement is often performed instead of repair due the technical difficulty of valve repair [18]. Valve repair typically requires use of cardiopulmonary bypass, a procedure which involves arresting the heart and emptying it of blood. The surgeon must try to imagine how the valve leaflets, and/or the chords that tether them, must be modified to make the valve close effectively after the heart is refilled with blood and pumping has been restored. Practice and experience are crucial for the development of the skills necessary to reliably repair mitral valves. Studies show that experienced surgeons at large clinical centers have a much better record of successful repairs, and valve replacement is often chosen instead of repair at low volume centers [8].
Surgical simulation has the potential to enable less experienced surgeons to effectively repair valves. This would allow many patients to undergo valve repair who would otherwise have undergone valve replacement. We propose that computer simulations of mitral valve closure can be used to help the surgeon plan effective repair strategies on a patient-specific basis. Under the proposed scheme, the geometry of a particular patient's valve would be extracted from medical images acquired prior to the date of surgery. The surgeon could then modify a computer model of that valve to reflect a particular surgical repair strategy and would use computer simulation to predict the closed state of the valve, indicating the effectiveness of that particular repair strategy. In this way, many potential repair strategies could be simulated and compared prior to the actual surgery, informing the surgeon to help choose the best strategy for a particular patient.
An important component of the proposed surgical simulation system is the computational model of the mitral valve. The proposed surgical simulation environment places two important requirements on the computational model. First, valve anatomy must be modeled in sufficient detailed to allow predictive modeling on a patient-specific basis. Second, the model must be able to compute the closed state relatively quickly. A surgeon may want to simulate ten or more surgical strategies for a given patient, so the time to simulate one valve closure must be on the order of minutes.
Several groups have developed finite element models of the mitral valve to study its function [6,11,23]. While these studies modeled important aspects of the complex behavior of the valve, their methods are not well-suited for the surgical simulation environment. They were based on averaged valve data, rather than subject-specific images, assuming symmetry of the leaflets through their midline and neglecting the branching structure of the chords. Another finite element study modeled the valve structures asymmetrically and obtained boundary conditions dynamically using implanted sonomicrometry crystals in an animal model [12]. Unfortunately their sonomicrometry method cannot be used clinically. All of these finite element models have execution times that are too slow for this surgical planning application.
In developing the computational model, several assumptions were made. First we assumed that a static loading state of peak systolic pressure was sufficient to assess valve competency in the model. A justification of this assumption is that the technique used at the end of surgery to test the repaired valve is to load the valve by injecting saline under static pressure [5]. This assumption allows us to ignore the complex interaction between blood flow and the valve structures during ventricular filling and ejection.
The second assumption concerns the role of tissue deformation in determining the shape of the closed valve. The valve leaflets undergo both conformational changes and deformation (tissue strain) in going from the open to the closed, loaded state. While the constitutive properties of valve leaflets and chords are known to be complex and to play a role in maintaining relatively low and uniform stress concentrations across the valve leaflets, we hypothesize that the conformational changes largely dictate whether the valve closes completely and that modeling the conformational changes along with a simple model of tissue properties will enable us to accurately predict the closed state given a particular valve geometry.
To meet these requirements, we have developed a computational model based on a mass-spring system, a method used in computer graphics to simulate the dynamics of fabric [16]. Mitral valve geometry is read directly from computed tomography (CT) data. This data is used to generate a triangular mesh. The mesh is treated as a system of masses connected by springs, and dynamics equations are used to evolve the closed state of the valve. The closed state predicted by the model is compared directly with images of the actual valve taken in the closed state.

Imaging
The mitral valves of two explanted porcine hearts were statically loaded with air via tubing inserted through the aorta, past the aortic valve, and into the left ventricle. The aorta was then cinched tightly around the tubing. To prevent air leakage through the coronary arteries, they were sutured closed. In order to supply air at low pressure with high accuracy, a circuit consisting of low-pressure regulators and electronic pressure sensors was constructed. The hearts were imaged in two different states using a micro-CT system (microCAT, Siemens, Munich, Germany): (1) with the mitral valve in the open position (no applied pressure), and (2) with the mitral valve in the closed position under typical porcine peak systolic pressure of 100 mmHg. Images were acquired at 100 μm isotropic voxel size. The volumetric CT image of the hearts were cropped to include only the mitral valve leaflets and chords. The resulting image of the valve was segmented, and an isosurface was fit to the data in Matlab (Mathworks, Natick, MA). The surface consists of an unstructured triangular mesh of points covering all surfaces of the leaflets and chords. The set of triangles comprising the atrial surface of the leaflets was isolated, and all chords that attach to either the free edge or the belly of the leaflets were approximated with line segments.

Model Structure
The dataset consisting of the triangulated mesh of the open valve leaflets along with the line segments representing the chords was used as the basis for a mass-spring model. All edges of triangles were treated as translational springs supporting only tension, and the mass of each triangular element (assuming finite thickness and known mass density) is treated as being lumped at the nodes. An example of a simple massspring mesh is shown in Figure 2. Spring constants for the springs comprising the valve leaflets were chosen using the following equation for approximating elastic membrane behavior with spring meshes [22]: (1) where k c is the spring constant for a given triangle side, E 2 is the two-dimensional Young's modulus for the leaflet tissue, the summation term represents the area of all triangles sharing side c, and the denominator is the squared length of side c. The two-dimensional Young's modulus is the product of Young's modulus and leaflet thickness. We assume uniform leaflet thickness of 1mm.

Figure 2
Example of a simple mass-spring mesh. All triangle sides are treated as translational springs, and mass is lumped at the nodes.
The stress-strain relationship for mitral valve leaflet tissue in known to be nonlinear, with a highly extensible pre-transitional region followed by a linear post-transitional region of much higher stiffness. See Figure 3. We approximated this relationship using a bilinear fit, with pre-and post-transitional stiffness of 100 and 6000 kPa and transition point of 25% strain [14]. Chord segments were also treated as springs supporting only tension, and spring constants were computed as 1-d Young's moduli based on chord length, cross-sectional area and Young's modulus for the chords [11]. Nodal mass was computed as the product of the nodal area (one third of the sum of the areas of triangles sharing that node), leaflet thickness and mass density.

Figure 3
Example of typical stress-strain curve observed in mitral valve leaflets. Young's modulus of the pre-transitional region, E pre , is the slope of the stress-strain curve at low strains, and Young's modulus of the post-transitional region, E post , is the slope at high strains. The transition point is denoted as ε*.

Model Dynamics
The dynamics of the mass-spring system can be expressed in state-space form as: (2) where x and v are vectors of nodal positions and velocities, respectively, M -1 is the inverse mass matrix (a diagonal matrix with the reciprocal of nodal mass on the main diagonal), and f is the vector of net nodal force due to springs and external forces. Implicit numerical integration is used because it allows larger integration step sizes and correspondingly faster simulations [1]. In order to use implicit integration, we discretized (2) using a second-order backward-difference formula as: where h is the integration time step. The net nodal force at step n+1 depends on the nodal positions at step n+1 making the set of equations nonlinear. It can be linearized by replacing f at step n+1 with a firstorder Taylor series approximation: Following a method used in a study simulating the behavior of cloth [3], (3) and (4) can be combined and expressed as the linear system: The Jacobian matrix expressing the partial derivative of the net force vector with respect to velocity is an N x N block matrix where N is the number of nodes in the system, and each block is 3 x 3, representing the three spatial coordinates. The forces due to springs as well as those due to applied pressure do not depend explicitly on nodal velocity, so their contributions are zero. Only the viscous damping term depends on nodal velocity, and its partial derivative yields -bI where b is the damping coefficient and I is the 3N x 3N identity matrix.
The Jacobian matrix expressing the partial derivative of the net force vector with respect to position is the same size as the Jacobian described above. In this case, the forces due to viscous damping and those due to applied pressure do not depend explicitly on position, so their contributions are zero. The forces due to the translational springs depend directly on nodal position, and their contribution to the Jacobian was evaluated analytically. For the translational spring between nodes i and j, elements of the Jacobian are computed as: (6) and (7) where (8) In this equation, l is the undeformed length of the spring between nodes i and j, {s x s y s z } T is the vector from node i to node j, and r = {s x s y s z }*{s x s y s z } T .

Solution Method
Equation (5) is a linear system where the first term on the left side is a sparse 3N x 3N matrix and the second term is a 3N x 1 vector of unknowns. All of the terms on the right side are 3N x 1 vectors which are known. It can be solved by inverting the sparse matrix. We used an iterative technique based on the method of conjugate gradients [2].
Points lying on the annulus as well as the locations where chords attach to the heart wall are treated as fixed (zero-displacement). However, both of these sets of points move considerably as the valve closesboth during physiological valve closure and during the passive loading that we use to image the closed valve. The closed shape of the valve leaflets is strongly dependent upon the locations of the annulus and chord attachment points, so it is important that we use their positions in the closed state for our simulations. To do so, we took a CT scan of the valve in the closed, loaded state then generated a mesh and identified the annulus and chord attachment points on the mesh. The annulus points were registered to those from the mesh of the open valve using the iterative closest point algorithm [24]. Points lying on the annulus of the open mesh were then linearly warped onto the annulus from the closed image, and the  points of attachments of the chords were moved directly to their positions measured in the image of the closed valve. All of the nodes in the mesh will be disturbed by the jump in positions of the annulus and chord attachments. To calculate their equilibrium state, the spring network was solved using a quasistatic approach. This was done by assembling the global stiffness matrix and solving it subject to zerodisplacement boundary conditions on the free edge of the leaflet and prescribed-displacement boundary conditions on the annulus and chord attachment points [13]. Then, during dynamic simulations, the constraint on the free edge of the leaflet is relaxed while annulus and chord attachment points are constrained to remain fixed.
Zero-displacement boundary conditions are implemented during simulations through use of the inversemass matrix appearing in (5). A particle i acted upon by springs but not subject to any displacement constraints will contribute the 3 x 3 diagonal matrix given by (1/m i )I to the main diagonal of the 3N x 3N inverse-mass matrix. However, we could prevent the velocity of the particle from changing by making the inverse-mass equal to zero, i.e., giving it an infinite mass. An infinite mass cannot be accelerated, so it effectively ignores all forces exerted on it. The zero displacement boundary conditions at the mitral valve annulus and at nodes where chords terminate in the heart wall are handled this way. Self-collisions of the leaflet were identified using a simple method based on proximity of vertices. Detected collisions were handled by inserting forces to render the collisions inelastic.

Model Parameters and Implementation
Some of the model parameters, such as constitutive properties of the tissues and applied transleaflet pressure, affect the closed shape of the valve at equilibrium. These parameters are assigned physically realistic values and are listed in Table 1. The remaining model parameters affect model dynamics and/or stability but not the closed shape of the valve, and those are assigned in order to minimize execution time and instability. The model was implemented in the Matlab programming language.

Results
Images from several stages of the simulation process for two different data sets are shown in Figure 4.  To quantitatively compare the closed state predicted by the model to the closed state generated from the image of the closed valve, the two surfaces were co-registered, again using the iterative closest point method based on vertices lying on the valve annuli. The error in the closed state predicted by the model is estimated by computing the magnitude of the distance between points on the closed image and their nearest points on the closed model. This distance is mapped to color and is plotted in Figure 5, with the error map on the left and right corresponding to the data sets in the top and bottom rows of Figure 4. The mean error across the surface was 1.7 mm for the error map on the left and 1.1 mm for the error map on the right. Maximum error was about 4 mm for both error maps.  The sensitivity of model results to changes in several important model parameters was evaluated. We define sensitivity, S, as: where Y is the measure of model accuracy, X is a parameter being tested, and X 0 is the value of that parameter used for our simulations and listed in Table 1. For Y, we use the mean error across the model surface. We approximate (9) as ΔY/ ΔX by increasing parameter X by 10%, repeating a simulation, and computing the resulting change in Y. Sensitivity to the constitutive properties of the leaflets is shown in Table 2.

Discussion
The goal of this study was to develop a simplified model of mitral valve mechanics specifically for use in surgical planning. There are three main requirements for the model. First, the model must represent the geometry of the valve structures in sufficient detail to allow patient-specific simulation. Second, the model must be able to simulate valve closure quickly and robustly. The third requirement concerns accuracy. Each of these requirements will be discussed below.
To produce models capable of conveying patient-specific anatomical detail, we produced dense meshes directly from images. Our imaging method provided high resolution and contrast and enabled us to acquire images under carefully controlled loading conditions. Micro-CT scans cannot be used to acquire images in the clinical setting because of the small bore diameter, and a patient's heart cannot be statically loaded for imaging. However, flat-panel volume CT can be used to image a human heart in vivo with similar resolution to our data [10]. Cardiac gating allows images to be captured at any point in the cardiac cycle, obviating the need for static loading.
Our mitral valve models were able to simulate one closing cycle in approximately 5 minutes, and significant speed gains can likely be made by implementing some of the bottleneck sections of the program in the C programming language. Further gains could be made by taking advantage of multiple CPU's or by using the GPU [19]. Simulations proved to be very robust. They were stable for all meshes that were tested, and stability was not affected by the quality of triangles in the mesh.
Parameter Sensitivity In choosing mass-spring modeling over finite element approaches, we have deliberately traded off some accuracy in the interest of speed and robustness. Finite element methods are based on continuum mechanics and can rigorously handle the anisotropy and nonlinearity that are known to characterize valve biomechanics [17]. Furthermore, they provide detailed analysis of stresses throughout the structures under load. However, our accuracy goals are more modest. At present, the surgeon hopes simply to create a mitral valve repair geometry that closes completely at peak load; our model is presented as a tool that could better inform surgeons as they try to understand the relationship between the geometry of the opened valve and its closed state. Analysis of stress concentrations throughout the leaflets is beyond the capability of the type of model presented here.
The model is able to predict many features of the closed state accurately and estimates the actual position of the closed leaflets with mean errors of 1.7 mm or less. By quantifying and plotting the error in the closed state predicted by the model, we can clearly see in which regions the model succeeds or fails to capture the actual behavior. For both mitral valves that we modeled, the maximum error of approximately 4 mm occurred in the middle of the leaflets. Two factors probably contribute to this error. First, by representing the mitral valve leaflets as isotropic, we neglect it strong orthotropic behavior, which is likely to play a role in determining leaflet shape. Second, for chords that attach to the free edge of a leaflet, we attach them at a single point on the edge, while, in reality, the chord inserts into the leaflet over a long overlapping region and imparts high stiffness to that leaflet in the direction of the chord.
The limited sensitivity analysis that we performed demonstrates that the accuracy of the model in predicting the closed state of the leaflets is not highly sensitive to the choice of leaflet properties. For example, a 1 kPa increase in the pre-transitional Young's modulus for the leaflets results in a decrease in model error of less than 1/100 th of a millimeter. It is desirable for our model to be relatively insensitive to leaflet properties because it indicates that we could have used any physiological values for leaflet properties (which are known to exhibit a large variance [14]) without significantly affecting our results.
It is important to note that closure of the valve leaflets is not the only metric of valve function, and hence quality of potential repair. One might also consider stress levels in the leaflets, a metric important for long-term durability of the valve. However, accurate simulation of the closed state is a good first-order criterion for valve function.

Conclusions
Our method of simulating closure of the mitral valve meets the requirements of surgical planning for valve repair. Simulations are fast and robust, and patient-specific models can be derived directly from images. Results are in reasonable agreement with images of the loaded valve. The relationship among the full set of model parameters need to be better understood, and the effect of changing the mesh density on speed and accuracy needs further investigation.