Patent Publication Number: US-9846765-B2

Title: System and method for patient specific modeling of liver tumor ablation

Description:
This application claims the benefit of U.S. Provisional Application No. 61/724,567, filed Nov. 9, 2012, the disclosure of which is herein incorporated by reference. 
    
    
     BACKGROUND OF THE INVENTION 
     The present invention relates to ablation planning, and more particularly to therapy planning and guidance based on a patient-specific model of liver tumor ablation using medical imaging data. 
     Ablation is one option for cancer treatment. In spite of recent advances in cancer therapy, treatment of primary and metastatic tumors of the abdomen remains a significant challenge. Hepatocellular carcinoma (HCC) for example is one the most common malignancies encountered throughout the world (e.g., &gt;1 million cases per year). In the U.S. alone, 1 in 153 individuals will develop HCC with reported 5-year survival rates of less than 15%. 
     For both primary liver cancer and hepatic metastases, liver resection (partial hepatectomy) is the current preferred option in patients having confined disease. In selected cases of early HCC, total hepatectomy with liver transplantation may also be considered. Unfortunately, less than 25% of patients with primary or secondary liver cancer are candidates for resection or transplantation, primarily due to tumor type, location, or underlying liver disease. Consequently, increasing interest has been focused on ablative approaches for the treatment of unresectable liver tumors. Rather than extirpation, this technique uses complete local in situ tumor destruction. A variety of methods have been employed to locally ablate tissue. Radiofrequency ablation (RFA) is the most commonly used, but other techniques are also used, including ethanol injection, cryo-therapy, irreversible electroporation, and microwave ablation. 
     The RFA procedure is performed by placing a probe within the target area, the tumor, in the liver parenchyma. Electrodes at the tip of the probe create heat, which is conducted into the surrounding tissue, causing coagulative necrosis at temperatures between 50° C. and 100° C. In addition to increasing the number of patients eligible for curative therapy of liver cancer in unresectable patients, local tissue ablation has a significant advantage as ablation may be performed using a minimally invasive approach, including percutaneously and laparoscopically. 
     To place the probe at the target location, the physician relies on intra-operative imaging techniques, such as ultrasound. However, the success of the procedure depends on the optimal placement of the probe and heat delivery. Different placements may have different results. The success of ablation is further challenged by the hepatic blood vessels that dissipate heat, thus potentially reducing RFA efficiency. 
     BRIEF SUMMARY OF THE INVENTION 
     The present invention provides a method and system for patient-specific modeling of liver tumor ablation. Embodiments of the present invention model the effects of ablation including heat diffusion, cellular necrosis, and blood flow through the vessels and liver. Starting from a preoperative medical image of a patient, such as a computed tomography (CT) or magnetic resonance imaging (MRI) image, the patient-specific geometry of the liver and venous system are automatically estimated. The vessel structure for the patient is accounted for as a heat sink in the model of biological heat transfer. The biological heat transfer model is coupled to a computational fluid dynamics solver to accurately take into account the effect of blood circulation on the dissipated heat. A cellular necrosis model is employed to simulate cell death due to over-heating and a simulated necrosis area can be visualized. 
     In one embodiment of the present invention, a patient-specific anatomical model of the liver and venous system of the liver is estimated from 3D medical image data of a patient. Blood flow in the liver and the venous system of the liver is simulated based on the patient-specific anatomical model. Heat diffusion due to ablation is simulated based on a virtual ablation probe position and its operating parameters and the simulated blood flow in the liver and the venous system of the liver. Cellular necrosis in the liver is simulated based on the simulated heat diffusion. A visualization of a simulated necrosis region and temperature map is generated. 
     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. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  illustrates a method for patient-specific modeling of liver tumor ablation according to an embodiment of the present invention; 
         FIG. 2  illustrates exemplary results for estimating a patient-specific anatomical model of the liver; 
         FIG. 3  illustrates an exemplary computational fluid dynamics (CFD) model of the hepatic venous circulatory system; 
         FIG. 4  illustrates an algorithm for implementing a computational model for simulating radio frequency ablation (ABA), according to an embodiment of the present invention; 
         FIG. 5  illustrates an exemplary visualization of a simulated necrosis area; 
         FIG. 6  illustrates an exemplary temperature map; 
         FIG. 7  illustrates analysis of spatial and temporal simulation results compared to an analytical solution; and 
         FIG. 8  is a high-level block diagram of a computer capable of implementing the present invention. 
     
    
    
     DETAILED DESCRIPTION 
     The present invention relates to patient-specific modeling and simulation of liver tumor ablation using medical imaging data. Embodiments of the present invention are described herein to give a visual understanding of the methods for patient-specific modeling and simulation using medical imaging data. 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 utilize a computational framework for patient-specific planning and guidance of radio-frequency ablation. Starting from a preoperative 3D medical image, such as a computed tomography (CT) image, the geometry of the liver and the venous system is automatically estimated using an efficient machine algorithm. A bio-heat equation is then solved to obtain a temperature profile throughout the liver over time. Tissue parameters are updated at every time step of the computation according to a cellular necrosis model. 
       FIG. 1  illustrates a method for patient-specific modeling of liver tumor ablation according to an embodiment of the present invention. The method of  FIG. 1  transforms medical image data representing a patient&#39;s liver anatomy to provide patient-specific simulation of liver tumor ablation. It is to be understood that the method of  FIG. 1  is not limited to the liver and can be similarly applied to other target organs as well. At step  102 , pre-operative 3D medical image data of at least of liver region of a patient is received. The pre-operative 3D image data can be acquired using any type of medical imaging modality, such as computed tomography (CT), three-dimensional rotational angiography, magnetic resonance imaging (MRI), ultrasound (US), positron emission tomography (PET), etc., provided that the liver is entirely visible in the medical image data. The medical image data can be received directly from an image acquisition device, such as a CT scanner, a C-arm image-acquisition device, an MRI scanner, or an US scanner, or the pre-operative cardiac image data can be received by loading previously stored cardiac image data of the patient. In a possible embodiment, a computed tomography, fluoroscopy, x-ray, or CT angiography system may be used to acquire the medical image data. The patient may ingest or be injected with a contrast agent that is generally opaque to x-rays. The contrast agent collects or is in the circulatory system, so vessels contrast with tissue. Any scanning sequence or approach may be used, such as CT angiography modes or rotational angiography to generate the pre-operative 3D medical image data. 
     At step  104 , an indication of a position of a virtual ablation probe and its operational parameters are received. In one embodiment, a user indication of a placement position of the ablation probe is received. The ablation is to occur near or in a lesion (e.g., tumor) or other region of tissue. Given anatomy, type of ablation, type of ablation device, or other limitations, the placement in the tissue may be limited. The user may indicate a possible placement of the ablation probe by selecting a location or locations in the pre-operative medical image data. For example, a location in a multi-planar reconstruction of the patient can be selected by a user using an input device, such as a mouse, touch screen, etc. 
     In a possible implementation, multiple positions of the ablation probe may be provided. For example, the user may indicate a sequence of placements for simulating sequential ablation operations or applications. In this case, the sequence is simulated by repeating steps  108 - 114  of  FIG. 1  for each sequential probe position using the results from the previous run. Using the modeling of cellular necrosis in step  112 , the changed tissue properties for various locations may be considered during the subsequent runs. Sequential placement may be used for larger tumors, where the single probe placement does not provide sufficient coverage of the thermal dose to the tumor. In another possible implementation, the user may indicate multiple placements for simulating ablation using multiple devices at a same time. The aggregated thermal dose is computed based on a single simulation run. 
     In an alternative embodiment, the user does not indicate placement. Instead, the position is selected automatically based on the image data, such as by identifying a center of the tumor. Various possible placements may be automatically identified and tested with separate simulations. 
     In addition to the position of the virtual ablation probe, the user may also input a spatial extent of the ablation probe, the type of ablation, duration, a desired dose, an indication of the spatial extent of the tumor, an indication of a location in the tumor, an amount of power for ablation, a type of ablation device, a sequence of power, and/or other characteristic of the ablation or tissue. Various inputs may be automated. Instead of user input, a processor provides the information. 
     At step  106 , a patient-specific anatomical model of the liver and the circulatory system in the liver is estimated from the 3D medical image data. The patient-specific anatomical model is a detailed anatomical model of the patient&#39;s liver parenchyma, tumors, hepatic vein, portal vein, and arteries. The 3D surface of the liver may be segmented automatically from the 3D medical image data. In one embodiment, the patient-specific anatomical model is semi-automatically estimated from the CT data using a graph theory approach, such as random walker segmentation. For each structure (parenchyma, tumors, hepatic vein, portal vein, arteries), the user defines seeds inside and outside the area of interest. Then, the random-walker algorithm automatically estimates the boundary of the structure. The process may be refined by the user, interactively, if necessary. Next, the resulting segmentations are merged into a multi-label mask image, which is then used to generate a tetrahedral multi-domain mesh. 
       FIG. 2  illustrates exemplary results for estimating a patient-specific anatomical model of the liver. As shown in  FIG. 2 , image  200  shows the liver segmentation  102  overlaid on a CT image. Image  210  shows the tetrahedral volume mesh generated from the liver segmentation, including the parenchyma  212 , hepatic vein  214 , portal vein  216 , and a segmented tumor  218 . 
     Returning to  FIG. 1 , at step  108 , blood flow in the liver is simulated based on the patient-specific anatomical model of the liver and venous system. The blood flow through the venous system of the liver acts as a heat sink in diffusing the heat applying by the ablation probe. The blood flow in the patient-specific anatomic model is simulated in order to provide a personalized modeling of the heat sink due to the blood flow. Locations of the vessels in the tissue region, the size of the vessels, and/or other vessel characteristics are used to model the heat sink characteristics in the region of interest. The characteristics may be extracted from the segmented vessel information in the patient-specific anatomical model of the liver. 
     The liver is treated as a porous media, so solving Darcy&#39;s law provides a velocity field throughout the organ that is used in the advection part of the heat transfer model. The blood velocity v inside the parenchyma is calculated according to Darcy&#39;s law: v=−k/(μe 2/3 )del·p, where p is the pressure within the parenchyma and μ is the dynamic viscosity of the blood flow. This amounts to solving the Laplace equation del·((−k/(μe 2/3 )del·p)=0. At the border of the liver, Neumann boundary conditions are employed as no or little flow leaks outside the liver, whereas at the tip of portal and hepatic veins, Dirichlet boundary conditions are applied. Since the pressures cannot be estimated in-vivo, a computational fluid dynamics (CFD) model of the hepatic venous circulation system is used to estimate the pressures.  FIG. 3  illustrates an exemplary CFD model of the hepatic venous circulatory system  300 . As shown in  FIG. 3 , arrows denote blood flow, circles denote portal vessel tips, and squares denote hepatic vessel tips. Let φ vcin  be the vena cava inflow, φ p  the portal vein inflow, and φ vc  be the vena cava outflow, which in virtue of conservation of mass is equal to φ vcin +φ p . The hepatic artery may be ignored or included as well. A small pressure is assumed to exist at the outlet of the vena cava (e.g., p 0 =1 mmHg). 
     The blood flow and pressure distribution are computed within the vena cava ( 302  in  FIG. 3 ) using three-dimensional computational fluid dynamics (e.g., using unsteady incompressible Navier-Stokes equations with viscous terms). The blood is modeled as a Newtonian fluid with pre-specified density (e.g., density=1.05 g/cm 3 ) and viscosity (e.g., viscosity=0.004 Pa s). The parameters may be different from one patient to another. A plug profile velocity field is applied at the inlets (squares in  FIG. 3 ), computed from the outflow φ p  and the cross-sectional area of each inlet. From the computational fluid dynamics computation, the downstream pressures p i   −  for each inlet of the vena cava is provided. The upstream pressure p +  at each of the portal vein tips (circles in  FIG. 3 ) is assumed constant such that the computed perfused flow through the vena cava inlets matches the three-dimensional computational fluid dynamics inlet flow profiles. Once p +  is estimated, the blood flow inside the portal vein ( 304  in  FIG. 3 ) is computed using the three-dimensional computational fluid dynamics solver. The effect of heat on the viscosity of the flow may be neglected. As a result, the flow related computation to simulate the blood flow in the liver can be performed before the simulation of heat diffusion. 
     Returning to  FIG. 1 , at step  110 , heat diffusion in the liver due to ablation is simulated based on the simulated blood flow and the location of the virtual ablation probe. The diffusion of heat over time in the liver tissue is simulated by calculating the diffusion or progression of temperature change over time. A temperature distribution field is solved as a function of time using a partial differential equation or other equation. The numerical solution of the bio-heat transfer equations yields an underlying temperature distribution field that varies in space and time. Computing heat diffusion in biological tissue amounts to solving the coupled bio-heat equations derived from the theory of porous media, where each elementary volume is assumed to comprise a tissue fraction and a blood fraction. Two main simplifications for solving the coupled bio-heat equations are the Pennes model and the Wulff-Klinger (WK) model. 
     In the Pennes model, the blood temperature is assumed constant, which holds close to large vessels, where blood velocity is high. The Pennes model may be expressed as:
 
(1−ε) p   t   c   t (δ T   t   /δt )=(1−ε) Q +(1−ε) del ·( d   t   del·T   t )+ H ( T   b0   −T   t ).  (1)
 
     In the Wulff-Klinger model, equilibrium between tissue and blood temperatures is assumed (Tt=Tb). This model is therefore well-suited small vessels, where blood velocity is low. The Wulff-Klinger model can be expressed as:
 
(1−ε) p   t   c   t (δ T   t   /δt )=(1−ε) Q +(1−ε) del ·( d   t   del·T   t )−ε p   b   c   b   v·del·T   t .  (2)
 
     In both equations, T is the temperature, Q is the source term and v is the blood velocity. Subscripts t and b refer to tissue and blood phases respectively. T b0  is the mean temperature, which is assumed to be constant, of the blood in the portal and hepatic veins. Definitions and exemplary values of the remaining model parameters are provided in Table 1, below. 
     
       
         
           
               
             
               
                 TABLE 1 
               
             
            
               
                   
               
               
                 Values of the parameters used in the simulation 
               
            
           
           
               
               
               
            
               
                 parameter 
                 description 
                 value 
               
               
                   
               
               
                 ρ b  · ρ t   
                 blood and tissue donation 
                 1.06 × 10 3  kg m −3   
               
               
                 c b   
                 blood heat capacity 
                 4.18 × 10 3  J(kg K) −1   
               
               
                 c t   
                 tissue heat capacity 
                 3.6 × 10 3  J(kg K) −1   
               
               
                 c t * 
                 tissue heat capacity in 
                 0.67 × 10 3  J(kg K) −1   
               
               
                   
                 dead cells 
                   
               
               
                 d b  · d t   
                 blood and tissue heat 
                 0.512 × (1 + 0.00161 × 
               
               
                   
                 conductivities 
                 (T − 310)) W(m K) −1   
               
               
                 H 
                 connective transfer coefficient 
                 24.4 × 10 4  W (m 3  K) −1   
               
               
                 ε 
                 blood volume fraction 
                 0.1 
               
               
                 κ 
                 permeability 
                 2.0 × 10 −2  m 2   
               
               
                 μ 
                 dynamic viscosity 
                 4.6 × (2.414 × 10 −5  × 10 247.8 / 
               
               
                   
                   
                 (T − 140) Pas   
               
               
                 k f   
                 forward rate constant 
                 3.33 × 10 −3  s −1   
               
               
                 k b   
                 backward rate constant 
                 7.77 × 10 −3  s −1   
               
               
                 T k   
                 parameter of cell state model 
                 40.5° C. 
               
               
                   
               
            
           
         
       
     
     The main difference between the Pennes and the Wulff-Klinger models lies in their cooling terms (i.e., the last term of the right-hand side). The former acts as a volumetric homogeneous heat sink whilst the latter accounts for the directional effect of blood flow on the tissue temperature field. Therefore, both equations can be easily implemented in a modular way to cope with tissue inhomogeneity, by using one or the other cooling term. 
     Current imaging techniques may not allow estimation of an accurate ratio between blood and liver tissue. Large vessels are clearly identified in the patient-specific anatomical model of the liver, but small capillaries are difficult to image. Therefore embodiments of the present invention solve the bio-heat equation by combining the Pennes and Wulff-Klinger models in a unified and modular framework to model large and small vessels. Assuming that blood vessels and the surrounding tissue are isolated from each other, the liver temperature T t  is solved using the diffusion equation p t c t (δT t /δt)=Q+del·(d t del·T t ) everywhere in the tetrahedral domain (i.e., everywhere in the domain of the patient-specific anatomical model). The cooling term H(Tb0−Tt)/(1−ε) is added when a tetrahedron belongs to a large vessel or artery (Pennes model) and the cooling term εp b c b v·del·T t  is added when the tetrahedron belongs to a small vessel or the parenchyma (Wulff-Klinger model). 
     The bio-heat and porous media models are solved using a finite element method (FEM). Test functions defined on linear tetrahedral elements, for instance, are used to discretize a weak form of the bio-heat equation (e.g., Garlekin&#39;s method). As a result, the heat diffusion is obtained by solving a linear system at every time step. Efficient implicit schemes may be used to have unconditional numerical stability, thus enabling large time steps for increased computational efficiency. A pre-requisite for FEM is the availability of detailed 3D tetrahedral, other meshes, or sampling of patient&#39;s liver. Once the anatomical model is generated, a virtual probe is placed and ablation is simulated using FEM. 
     In one example implementation of the finite element method on a tetrahedral mesh, the discretized equation: M +KU=B is solved, with U the temperature, M the mass matrix, K the stiffness matrix for the diffusion, the reaction and the advection part, and B the vector including boundary conditions. Diffusion, reaction and advection terms are implemented in a modular way, but may be combined. Neumann boundary conditions are used at the organ boundary. The heat source term can be modeled through a Dirichlet boundary condition of 100° C. at the location of the virtual probe. The blood going through the liver is treated as a constant temperature of 37° C. For the advection term, a three-dimensional characteristic Galerkin scheme is used for numerical stability. First-order implicit Euler time discretization is employed. The resulting matrix system Ax=B is solved using the conjugate gradient iterative algorithm. For the computational fluid dynamics, a full 3D Navier-Stokes viscous solver, expressed in an Eulerian framework which embeds the domain boundary using a level set representation of the segmented vessel is used. 
     Returning to  FIG. 1 , at step  112 , cellular necrosis in the liver is simulated based on the simulated heat diffusion. Heat from an ablation probe can cause cellular necrosis in the liver. In one embodiment, the tissue damage is modeled as using a three state model. The model computes the variation of concentration of live cells (A), vulnerable cells (C), and dead cells (D) over time according to the state equation: 
               A   ⁢           ⁢     •     k   b       k   f       ⁢           ⁢   V     ⁢     →     k   f       ⁢   D         
where k f  and k b  are the rates of cell damage and recovery, respectively. In particular, k f  depends on the temperature through the equation: k f =k f   − e T/T     k    (1 A) with k f   −  being a scaling constant and T k  a parameter that sets the rate of the exponential increase with temperature. This equation results in three simultaneous ordinary differential equations, solved at each node of the computational domain, yielding spatially-varying cell state field. Any temporal integration scheme may be used (e.g., first order Euler explicit or first order Euler implicit). In an exemplary implementation, the initial conditions can be chosen as A=0.99, V=0.01 and D=0.00. Exemplary values for the parameters, k b , k f   − , and T k  are shown in Table 1, above. Additionally, dead or damaged cells do not have the same heat capacity ct as live cells. This phenomenon can affect heat transfer and ablation extent. Accordingly, at each time step of the computation, the stated of the cell is updated at each point in the patient-specific anatomical model of the liver, given the cell&#39;s current temperature. The corresponding heat capacity ct is used if healthy or vulnerable, dead cells have a fixed heat capacity of c t *.
 
     Other cell death models may be used as well. For example, in another embodiment, the tissue damage can be modeled using a survival fraction index and an Arrhenius-like equation. The tissue damage model is experimentally determined, extracted from studies, or created from theory. The survival fraction index indicates the amount of live versus dead cells in a region as a function of temperature. The survival fraction index is, calculated from the current temperature using an Arrhenius-like model. More precisely, the survival fraction index is the ratio of viable cells to the total number of cells in a region of interest. When the survival fraction index is &lt;&lt;1, cells in this spatial region are dead. The Arrhenius equation models the relationship of the survival fraction with respect to temperature. 
       FIG. 4  illustrates an algorithm for implementing a computational model for simulating radio frequency ablation (ABA), according to an embodiment of the present invention. The algorithm of  FIG. 4  can be used to implement steps  106 - 112  of  FIG. 1 . As illustrated in  FIG. 4 , at  402 , the patient-specific model of liver anatomy is estimated. At  404 , compute 3D blood flow φi and inlet pressures p in hepatic vein are computed. At  406 , portal vein upstream pressure p+ is computed by fitting Darcy&#39;s model to the hepatic vein downstream pressures p and flow φi. At  408 , 3D blood flow in portal vein is computed. The blood flow and pressures in the hepatic vein and the portal vein provide Dirichlet boundary conditions for computing the blood velocity inside the parenchyma. Once the blood flow is simulated in the whole patient-specific model of the liver, operations  410  and  412  are performed at each time step (t) between while t&lt;t end . At  410 , the temperature T using Pennes model in the large vessels and Wulff-Klinger model elsewhere (e.g., parenchyma). At  412 , the cell-state is updated by using the cellular necrosis model based on the simulated temperature. 
     Returning to  FIG. 1 , at step  114  temperature maps and a visualization of a simulated necrosis area are output. For example, the temperature maps and simulated necrosis are visualization can be output by displaying the temperature maps and simulated necrosis are visualization on a display of a computer system. 
     The proposed framework has been evaluated on one patient-specific geometry extracted from pre-operative CT images. The virtual probe is placed in within the middle of the tumor. Cells within a sphere around the probe tip are heated at 100° C. to mimic the topology of the probe. The protocol to be clinically used is simulated as: heating for 10 minutes, stopping the heating for the next 5 minutes, and repeat the process by increasing the sphere diameter from 10 mm to 20 mm and finally 30 mm. Using the Pennes Model, the model is able to predict the cooling effect of the big veins. The distribution of the temperature over the liver and a simulated necrosis area are returned.  FIG. 5  illustrates an exemplary visualization of a simulated necrosis area. As shown in  FIG. 5 , the simulated necrosis area  502  is visualized by visualizing all portions of the tetrahedral mesh having a dead cell state. Alternatively, necrotic areas can be visualized using iso-contours of the survival fraction index. A ground truth cellular necrosis area  504  from a post-operative MRI is overlaid for comparison. As can be seen in  FIG. 5 , the predicted necrotic extent from the simulated necrosis area  502  qualitatively agrees with the ground truth cellular necrosis  504  that is actually observed in the patient. 
     The temperature maps show the spatial and temporal variation of the temperature in the liver over time. In one embodiment, the liver temperature can be visualized using color maps and isothermal contours that indicate the heat of particular regions in the liver.  FIG. 6  illustrates an exemplary temperature map. As shown in  FIG. 6 , a temperature map  600  showing temperature at a time step in the simulation is mapped to CT cross-section images. 
     The present inventors evaluated the computational model for simulating radiofrequency ablation in the liver by comparing its behavior on a regular tetrahedral beam mesh where the 3D analytical solution of an instantaneous point source of mass M released at location p 0 =(x 0 ,y 0 ,z 0 ) and time t 0  with v=(u,v,w) f the advection-diffusion equation ∂T/∂T+v·∇T=∇ (a∇ T) is known: 
               T   ⁡     (     p   ,   t     )       =       M       [     4   ⁢     π   ⁡     (     t   -     t   0       )       ⁢   d     ]       3   /   2         ⁢       exp   ⁡     [     -              p   -     p   0     -       (     t   -     t   0       )     ⁢   v            2       4   ⁢           ⁢     d   ⁡     (     t   -     t   0       )             ]       .             
The following values were used: M(° C./mm 3 )=100, d(mm 2 /s)=100, t 0 =0 s, and v(mm/s)=(0.001,0,0) (values of the same order of magnitude as for the patient simulation). The source term is set at the center of the mesh to minimize the effects of the boundary conditions. At time t=0 s, the temperature values were initialized at each vertex with the analytical solution at time t=0.05 s, and the temperature was analyzed at several points of the mesh.  FIG. 7  illustrates analysis of spatial and temporal convergence of simulation results compared to the analytical solution. Image  700  shows spatial convergence analysis comparing temperature values of the analytical solution  702  over time with temperature values calculated using simulations with mesh resolutions of 2 mm ( 704 ), 5 mm ( 706 ), and 10 mm ( 708 ). Image  710  shows temporal convergence analysis comparing temperature values of the analytical solution  712  over time with temperature values calculated using simulations with time steps of 0.005 s ( 714 ), 0.05 s ( 716 ), and 0.5 s ( 718 ). As expected, the more accurate smaller the time step and spatial resolution of the mesh, the more accurate the simulated solution is compared to the analytical solution. As shown in  FIG. 7 , a time step of 0.5 s and a resolution of 5 mm provides a good compromise between accuracy and computational cost.
 
     Embodiments of the present invention provide a patient-specific model for simulating liver tumor ablation. Embodiments of the present invention simulate heat propagation and cellular necrosis based on patient-specific anatomic model estimated from medical image data of the patient, and takes into account the heat sink effect of blood vessels and porous circulation in the liver. Embodiments of the present invention provide visualization of temperature distribution, and thus lesion propagation, as well as cellular necrosis, allowing for improved planning and guidance for liver tumor ablation procedures. 
     In a possible implementation, an optimal probe location may be automatically determined using the above described simulation framework. Rather than rely on user feed back, the temperature diffusion and corresponding necrosis associated with different locations and/or other variance in ablation can be compared to the region to be treated. For example, the position, power level, and sequence of application resulting in the greatest, sufficient, or complete coverage of the tumor and least necrosis of healthy tissue may be found by simulating different combinations. By successively computing the thermal dose for a number of needle placements, the placement with maximum coverage on the tumor can be found. In a successive optimization loop to solve for the probe location, orientation, and/or other characteristic, the combination for maximum tumor thermal dose coverage is found. The thermal dose coverage may be measured as the area defined by temperatures higher than a threshold (e.g., 70 degrees Celsius). Other measures may be used, such as a weighted measure that most strongly weights necrosis of the tumor while considering avoidance of necrosis outside the tumor. 
     The above-described methods for patient-specific modeling and simulation of liver tumor ablation 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  FIG. 8 . Computer  802  contains a processor  804 , which controls the overall operation of the computer  802  by executing computer program instructions which define such operation. The computer program instructions may be stored in a storage device  812  (e.g., magnetic disk) and loaded into memory  810  when execution of the computer program instructions is desired. Thus, the steps of the methods of  FIGS. 1 and 4  may be defined by the computer program instructions stored in the memory  810  and/or storage  812  and controlled by the processor  804  executing the computer program instructions. An image acquisition device  820 , such as a CT scanning device, C-arm image acquisition device, MR scanning device, Ultrasound device, etc., can be connected to the computer  802  to input image data to the computer  802 . It is possible to implement the image acquisition device  820  and the computer  802  as one device. It is also possible that the image acquisition device  820  and the computer  802  communicate wirelessly through a network. The computer  802  also includes one or more network interfaces  806  for communicating with other devices via a network. The computer  802  also includes other input/output devices  808  that enable user interaction with the computer  802  (e.g., display, keyboard, mouse, speakers, buttons, etc.). Such input/output devices  808  may be used in conjunction with a set of computer programs as an annotation tool to annotate volumes received from the image acquisition device  820 . One skilled in the art will recognize that an implementation of an actual computer could contain other components as well, and that  FIG. 8  is a high level representation of some of the components of such a computer for illustrative purposes. 
     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.