The present invention relates to the computation of advanced cardiac parameters from medical data and images by using a patient-specific multi-physics fluid-solid models of the heart.
The heart is the pump of life. The periodic contraction and relaxation cycle of the heart ensures blood circulation throughout the body. Heart diseases can affect the ventricle shape, the integrity of the muscle, the electrophysiology of the heart, etc., but the end result is incorrect or inadequate blood circulation in a patient. Hence, the ultimate goal of any cardiac therapy is to correct the blood flow, by acting on the various aspects of the cardiac system that may be affected, directly or indirectly, by existing pathologies in the patient.
Cardiac disease diagnosis and therapy planning are made difficult by the large variability in heart diseases. Every patient is unique. The variability observed among patients with the same disease is such that population-wise guidelines can be sub-optimal in terms of diagnosis, therapy outcome and complications for a specific patient. For example, patients with acute myocardium infarction can present with variable scar morphology and extent, which can influence the outcome of various cardiac therapies. Accordingly, a framework to assess the current state of the heart and the optimal therapy for a specific patient is desirable. Such a framework should be integrative and comprehensive, considering all major aspects of heart function, and focus on the end-result of any cardiac therapy, which is the restoration of normal blood flow.
The present invention provides a method and system for computing advanced cardiac measurements from clinical images and data using a patient-specific multi-physics model. Embodiments of the present invention provide a simulation framework that can be utilized to estimate advanced cardiac measurements not accessible through standard clinical devices and optimize various cardiac therapies, such as open-heart therapies or electrophysiology therapies, in terms of cardiac blood flow. Embodiments of the present invention utilize a fluid-structure interaction framework that couples a multi-scale model of cardiac electrophysiology, a multi-scale model of cardiac biomechanics, and a computational fluid dynamics based model of cardiac blood flow together in order to estimate intrinsic cardiac parameters, such as tissue stiffness, active stress, electrical diffusion, etc., and predict heart blood flow after a cardiac therapy. Medical imaging data and non-imaging clinical data of a patient are used to adjust free parameters of the simulation framework, such that the patient-specific heart function is accurately reproduced in the computer simulation. These parameters are then used as advanced cardiac measurements by the user. Therapies can then be simulated by altering the patient-specific model in terms of heart shape to simulate surgical therapies, or dynamics to simulate electrophysiology therapies.
In one embodiment of the present invention, a patient-specific anatomical model of left and right ventricles is generated from medical image data of a patient. A patient-specific computational heart model is generated based on the patient-specific anatomical model of the left and right ventricles and patient-specific clinical data. The computational model comprises an electrophysiology component, a biomechanics component and a hemodynamics component. The patient-specific computational heart model is generated by determining initial patient-specific parameters of an electrophysiology model, determining initial patient-specific parameters of a biomechanics model, determining initial patient-specific computational fluid dynamics (CFD) boundary conditions, performing a coupled fluid-structure interaction (FSI) simulation using the initial patient-specific parameters of the electrophysiology model, the initial patient-specific parameters of the biomechanics model, and the initial patient-specific CFD boundary conditions, and refining the initial patient-specific parameters of the electrophysiology model, the initial patient-specific parameters of the biomechanics model, and the initial patient-specific CFD boundary conditions based on the coupled FSI simulation. The refined parameters define a set of advanced cardiac measurements. At least one cardiac therapy is simulated using the patient-specific computational heart model.
These and other advantages of the invention will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings.
The present invention relates to a method and system to compute advanced cardiac measurements and plan cardiac therapies using a patient-specific multi-physics fluid-solid heart model that is personalized for a patient based on preoperative medical image data, such as MRI and/or Ultrasound data, and non-imaging clinical data, such as ECG and pressure measurements. Embodiments of the present invention are described herein to give a visual understanding of the cardiac therapy simulation methods. A digital image is often composed of digital representations of one or more objects (or shapes). The digital representation of an object is often described herein in terms of identifying and manipulating the objects. Such manipulations are virtual manipulations accomplished in the memory or other circuitry/hardware of a computer system. Accordingly, is to be understood that embodiments of the present invention may be performed within a computer system using data stored within the computer system.
Embodiments of the present invention provide a simulation framework that can be utilized to estimate advanced measurements, which correspond to the patient-specific values of the simulation framework, and to optimize various cardiac therapies, such as open-heart therapies or electrophysiology therapies, in terms of cardiac blood flow. Embodiments of the present invention utilize a fluid-structure interaction framework that couples a multi-scale model of cardiac electrophysiology, a multi-scale model of cardiac biomechanics, and a computational fluid dynamics based model of cardiac blood flow together in order to predict heart blood flow after a cardiac therapy. Medical imaging data and non-imaging clinical data of a patient are used to adjust free parameters of the simulation framework, such that the patient-specific heart function is accurately reproduced in the computer simulation. The estimated parameters constitute an output of the system, and are referred to herein as “advanced cardiac measurements”. Therapies can then be simulated by altering the patient-specific model in terms of heart shape to simulate surgical therapies, or dynamics to simulate electrophysiology therapies.
For each of the LV and the RV, the heart chamber segmentation can be formulated as a two-step learning problem: anatomical structure localization and boundary delineation. In an advantageous embodiment, marginal space learning (MSL) can be used to apply machine learning to 3D object detection. The idea of MSL is not to learn a monolithic classifier directly in the full similarity transformation parameter space but to incrementally learn classifiers on marginal spaces. In particular, the detection of each heart chamber can be split into three problems: position estimation, position-orientation estimation, and position-orientation-scale estimation. A separate classifier is trained based on annotated training data for each of these estimation problems. The classifiers in the lower dimensional marginal spaces are used to prune the searching space efficiently. This object localization stage results in an estimated transformation (position, orientation, and scale) of the object (e.g., heart chamber).
After automatic object localization, the mean shape model of the object is aligned with the estimated transformation to get a rough estimate of the object shape. The shape is then deformed locally to fit the object boundary. Active shape models (ASM) can be used to deform an initial estimate of a non-rigid shape under the guidance of the image evidence and the shape prior. However, a non-learning based generic boundary detector, as used in conventional ASM applications, does not work effectively in heart chamber deformation due to the complex background and weak edges. Instead, a learning based boundary detector can be used to exploit more image evidences to achieve a robust boundary detection. Additional details regarding MSL-based heart chamber segmentation are described in U.S. Pat. No. 7,916,919, issued Mar. 29, 2011, and entitled “System and Method for Segmenting Chambers of a Heart in a Three Dimensional Image”, United States Published Patent Application No. 2010/0040272, and United States Published Patent Application No. 2012/0022843, which are incorporated herein by reference.
At step 304, the patient-specific LV and RV models are fused into a single bi-ventricular myocardium volumetric mesh. In a possible implementation, the LV and RV anatomies extracted can be fused together. The resulting closed surface is used to create a volumetric, tetrahedral mesh on which vertices are tagged into surface zones according to the underlying anatomy.
At step 306, spatial information is mapped onto the bi-ventricular myocardium mesh. Spatial information, such as scars, grey zones, and fibrosis can be identified in images, such as late delayed-enhancement MR images and mapped onto the bi-ventricular myocardium mesh. For example, scar locations and extent can be segmented in delayed-enhancement MR images using the method described in Dikici et al, “Quantification of Delayed Enhancement MR Images”, In Proc. MICCAI 2004, LNCS 3216, pp 250-257, 2004, which is incorporated herein by reference. The scar information is mapped onto the bi-ventricular myocardium mesh by tagging the tetrahedral elements that lie within the segmented scar regions. This spatial information is important to simulate the electrical wave around scars, in particular for wave-reentry assessment, but also the impaired contractility due to dead tissue.
At step 308, model of fiber orientation is generated on the bi-ventricular myocardium mesh. In one embodiment, in-vivo diffusion tensor (DT) MR images of the patient's cardiac fibers are directly mapped to the anatomical model through image registration. In this case, the DT MR image is non-linearly registered to the medical image in which the LV and RV models are detected. The resulting transformation is used to deform the tensor field in the DT MR image towards the anatomical model. The Finite Strain method, the details of which are described in Peyrat et al., “A Computational Framework for the Statistical Analysis of Cardiac Diffusion Tensors: Application to a Small Database of Canine Hearts”, IEEE TMI, 26(11):1500-1514, 2007, which is incorporated herein by reference, is used to reorient the tensors once the tensors are registered to the anatomical model.
In another embodiment, if no in-vivo DT MRI are available, the model of fiber orientation may be computed directly from the anatomical model using a rule-based method. A generic model of myocardium fiber architecture that includes fiber and fiber sheets is computed. A rule-based strategy is followed to generate the fiber architecture to cover the entire bi-ventricular myocardium from apex to valves. Below the basal plane, which is identified automatically using point correspondences of the initial triangulations of the anatomical model, the fiber elevation angle α, i.e. their angle with respect to the short axis plane, varies linearly across the myocardium, from −70 on the epicardium to +70 on the endocardium. Similarly, the sheet direction, which is defined by the angle β with respect to the outward transmural axis, varies transmurally from +45 on the epicardium to −45 on the endocardium. α and β are computed for each point of the volumetric bi-ventricular myocardium mesh between the apex and basal plane based on the geodesic distance to the endocardia and epicardia identified by the facet tags: α=(depiαendo+dendoαepi)/(dendo+depi), where depi, dendo, αepi, and αendo are the distances and angles at the endocardium and epicardium, respectively. The fiber and sheet orientations are then fixed around each valve. In particular, fibers are longitudinal around the aortic valve and tangential around the mitral, tricuspid, and pulmonary valves, and sheet normals are oriented towards the barycenter of the valves. The local orthonormal basis is then interpolated from the basal plane to the valve, first by following the myocardium surface, then throughout the myocardium thickness. For orthonormality preservation, the interpolation can be performed using a Log-Euclidean framework.
Once the patient-specific anatomical model is generated, then patient-specific anatomical model is provided as an input to the computational heart model solver module (step 220 of
Referring to
where v(t) is the trans-membrane potential, h(t) is a gating function, K is the anisotropic diffusivity matrix, Jstim(t) is an external stimulus (a CRT lead for instance), and Tin, Tout, Tclose, Topen are free parameters that are directly related to the shape of the action potential. The model is paced at the septum with a user-defined frequency, which can be estimated from ECG signals.
In an advantageous implementation, the electrophysiology (EP) model can be extended to consider Purkinje fibers. In general, the anatomical structures of the Purkinje network cannot be extracted from medical imaging modalities. In order to mimic the fast propagation of electrical waves along the Purkinje fibers while keeping the problem tractable, embodiments of the present invention implement a “macro-scale” electrophysiology model of Purkinje fibers. For this model, it is assumed that the Purkinje fibers are evenly distributed on the entire endocardial surface. Thus, the Purkinje fibers can be modeled as a “continuous” membrane rather than a complex network. The above-described Mitchell-Schaeffer equations can be solved on this membrane using two-dimensional finite element formulation. In order to take into account the complex three dimensional shape of the endocardial surface, the finite element formulation is first built on each local, planar coordinate and then is rotated back into the global coordinate system to perform the final assembly. The fast propagation is implemented by setting a high diffusion coefficient KPurkin. The coupling between the “macro-scale” Purkinje fiber model and the myocardium is implemented through boundary conditions. Here, one way coupling is employed to link the Purkinje fiber with the myocardium. At each time step, the electrophysiological equation (Equation 1) is first solved on the endocardial triangulation defining the domain of the macro-scale Purkinje fiber. Since the endocardial triangulation has one-to-one mapping with the tetrahedral vertices at the endocardium, Dirichlet boundary conditions can be specified for the myocardium so that the potential at these tetrahedral nodes are set equal to the corresponding triangulation nodes. The equations are then solved on the entire myocardium. Additional details regarding the electrophysiology model are described in U.S. patent application Ser. No. 13/754,174, entitled “Method and System for Patient Specific Planning of Cardiac Therapies on Preoperative Clinical Data and Medical Images”, filed Jan. 30, 2013, which is incorporated herein by reference.
In order to personalize the electrophysiology model to a specific patient, the parameters of the electrophysiology model are adjusted using inverse problem methods based on the observed clinical data and medical image data for the patient. ECG and endocardial mapping are used to adjust tissue diffusivity K and action potential duration, mainly controlled by the parameter τclose, at every node of the anatomical model. For example, the simulated QRS (heartbeat) can be aligned with measured QRS by adjusting tissue diffusivity K. If endocardial mappings are available, inverse problem methods can be used to estimate tissue conductivity and τclose by comparing the simulated isochrones and potentials with measurements. Additionally, dynamic images can be used to adjust electrophysiology overall synchrony by minimizing the difference of the cardiac motion computed from the electrophysiology model and the observations, and bundle branch blocks can be simulated by modifying the activation times of the left and right ventricles.
For computational efficiency, a marginal approach is employed to estimate the initial EP parameters.
Returning to
MÜ+C{dot over (U)}+KU=F
c
+F
p
+F
b. (2)
In Equation 2, U(t) is the position of the node of the anatomical model, M is a diagonal mass matrix (mass density ρ=1.070 g/mL), and C is a Rayleigh damping matrix. K is the stiffness matrix of myocardium tissue. In a possible embodiment, the stiffness matrix K can be modeled using co-rotational transverse isotropic elasticity to increase computational efficiency while coping with large deformations and tissue anisotropy. The model is controlled by the Young Modulus E and the Poisson ratio v, which is set to 0.488 to model tissue incompressibility. In other possible embodiments, various types of non-linear hyperelastic models, such as a Mooney-Rivlin model or a Holzapfel-Ogden model, can be used to model the passive tissue properties of the myocardium tissue. Additional details regarding modeling the passive tissue properties are described in U.S. patent application Ser. No. 13/754,174, entitled “Method and System for Patient Specific Planning of Cardiac Therapies on Preoperative Clinical Data and Medical Images”, filed Jan. 30, 2013, which is incorporated herein by reference.
Depolarized cardiac cells contract according to an active force Fc(t) controlled by the trans-membrane potential v(t). That force increases exponentially as the cell depolarizes, at a rate +kATP related to the rate of ATP binding, to reach a maximum contraction strength σ0 that is controlled by the user. When the cell repolarizes, the force fades away exponentially at a rate −kRS related to the rate of ATP release. Additional details regarding modeling the active force Fc(t) of cell contraction are described in U.S. patent application Ser. No. 13/754,174, entitled “Method and System for Patient Specific Planning of Cardiac Therapies on Preoperative Clinical Data and Medical Images”, filed Jan. 30, 2013, which is incorporated herein by reference.
The base of the anatomical model is attached using stiff springs Fb(t) to mimic the compliance of the atria and arteries. A pericardium constraint can also be applied to keep the heart fixed in space. Nodes of the anatomical models that go outside the allowed region of the pericardium bag are moved back by applying a spring force, also gathered in the global vector Fb(t). Additional details regarding the base stiffness parameter and the pericardium constraint are described in U.S. patent application Ser. No. 13/754,174, entitled “Method and System for Patient Specific Planning of Cardiac Therapies on Preoperative Clinical Data and Medical Images”, filed Jan. 30, 2013, which is incorporated herein by reference. Fp(t) is a pressure applied on the endocardial surface of each ventricle resulting from the blood flow through the ventricles. That pressure is directly computed using the CFD solver described below.
In order to personalize the biomechanical model, the parameters of the biomechanical model (E, σ0) are adjusted from the dynamics images by comparing the observed cardiac dynamics with the simulated motion of the ventricles. For example, the active contraction strength σ0 can be adjusted from the motion observed during systole, and the tissue stiffness E can be adjusted to match the dynamics observed during passive systole.
Returning to
According to an advantageous embodiment of the present invention, the following boundary conditions are used, in order to achieve personalized CFD computations:
WALLS—The ghost fluid method is used to impose Dirichlet boundary conditions for the velocity and a finite volume cell flux balance using Neumann boundary conditions for the pressure;
INLETS—The mitral and tricuspid valves are used to impose specialized boundary conditions during all cardiac phases, as follows: during systole and isovolumetric phases, they are closed and the boundary conditions are set to no-slip wall; during diastole they are open, but their closed geometric position is used to impose Dirichlet boundary conditions for the personalized atrial pressures;
OUTLETS—The aortic and pulmonary valves are used to impose specialized boundary conditions during all cardiac phases, as follows: during diastole and isovolumetric phases, they are closed and the boundary conditions are set to no-slip wall; during systole they are open, but their closed geometric position is used to impose Dirichlet boundary conditions for the personalized aortic pressures. Such pressures can be obtained either from patient-specific clinical measurements, or from generative boundary conditions, based on a lumped model of arterial pressure like Windkessel models. Parameters of the Windkessel model can be personalized from pressure catheterization and arterial flow.
Valve kinematics can be modeled using a geometric approach as they are often difficult to visualize from patient-specific data. Each valve is defined as a triangulated patch prescribed by an annulus curve automatically detected on the outflow and inflow tracts of the ventricles. Machine learning approaches can be used to estimate the annuli from the images. Owing to point correspondence, the annuli can be tracked over-time, to determine the dynamics of the valve. Valve closure and opening are finally modeled by a ring whose inner diameter varies over time. Opening and closure timing are synchronized with patient data.
CFD boundary conditions are personalized from clinical data directly. Measured inflows and outflows are provided as input to the model. This data is useful in recovering patient-specific aortic/pulmonary pressures, which is done at every FSI simulation step by enforcing that the simulated flow due to endocardial pressurization matches the measured flow.
Returning to
Returning to
Returning to
Returning to
Pressure catheterization was available for this patient. Hence, the free parameters of the 3-element Windkessel model of the aorta were adjusted using an exhaustive search approach. Atrial pressures were set constant in this first experiment. Fluid-related parameters, such as flow and pressure are lumped at this stage. As no regurgitations were visible, arterial flow was equal to the opposite of the volume variation of the ventricle.
Although above described methods are described as being implemented with two-way coupling between the electro-mechanical model and the CFD simulation, the present invention is not limited thereto, and it is possible to implement one-way coupling as well. In the case of one-way coupling between the electro-mechanical model and the CFD simulation, myocardium function can first be simulated entirely using only electro-mechanics models, with blood flow being modeled with lumped parameters. The resulting myocardium displacements and velocities are then given to the CFD solver.
The above-described methods for generating a patient-specific anatomic model of the heart, generating a patient-specific computational model of the heart, and performing patient-specific cardiac therapy simulations can be implemented on a computer using well-known computer processors, memory units, storage devices, computer software, and other components. A high-level block diagram of such a computer is illustrated in
The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.
This application claims the benefit of U.S. Provisional Application No. 61/593,500, filed Feb. 1, 2012, the disclosure of which is herein incorporated by reference.
Number | Date | Country | |
---|---|---|---|
61593500 | Feb 2012 | US |