|Home | About | Journals | Submit | Contact Us | Français|
Elevating the temperature of cancerous cells is known to increase their susceptibility to subsequent radiation or chemotherapy treatments, and in the case in which a tumor exists as a well-defined region, higher intensity heat sources may be used to ablate the tissue. These facts are the basis for hyperthermia based cancer treatments. Of the many available modalities for delivering the heat source, the application of a laser heat source under the guidance of real-time treatment data has the potential to provide unprecedented control over the outcome of the treatment process [7, 18]. The goals of this work are to provide a precise mathematical framework for the real-time finite element solution of the problems of calibration, optimal heat source control, and goal-oriented error estimation applied to the equations of bioheat transfer and demonstrate that current finite element technology, parallel computer architecture, data transfer infrastructure, and thermal imaging modalities are capable of inducing a precise computer controlled temperature field within the biological domain.
Today, cancer is the second most common cause of death in the United States. In 2006, about 564,830 Americans are expected to die of cancer. Among the most common forms of the disease are cancer of the prostate and breast, accounting for 10% and 15% of all cancer related deaths in men and women, respectively. Current detection technology has led to a 5-year relative survival rate of greater than 90% for prostate and breast cancer . Given these survival rates, the quality of life of the patient post-treatment is a very important factor. Traditional medical treatment of early stage prostate and breast cancer, including radical prostatectomy, pelvic lymphadenectomy, lumpectomy, and mastectomy are major surgical procedures and are associated with many surgical complications. Ablation therapies delivered under various treatment modalities are a very promising, minimally invasive, alternative to standard treatment and show significant potential as an effective cancer treatment to eradicate the disease while maintaining functionality of infected organs and minimizing complications and relapse.
The physical basis for thermal therapies is that exposing cells to temperatures outside their natural environment for certain periods of time can damage and even destroy the cells. However, one of the limiting factors in all forms of ablation therapy, including cryotherapy, microwave, radio-frequency, and laser, is the ability to control the energy deposition to prevent damage to adjacent healthy tissue .
Recent advances in imaging technology allow the imaging of the geometry of tissue and an overlaying temperature field using MRI and new MRTI technology. MRTI has the desirable properties of excellent soft-tissue contrast, the ability to provide fast, quantitative temperature imaging in a variety of tissues, and the capability of providing biologically relevant information regarding the extent of injury immediately following therapy . Image guidance [18, 20] has the potential to facilitate unprecedented control over bioheat transfer by providing real time treatment monitoring through temperature feedback during treatment delivery. A similar idea using ultrasound guided cryotherapy has been studied and shows good results . However, the hypothermic regime requires temperature differentials of 70° from normothermia to ablate the tissue. Tissue ablation in the hyperthermic regime requires substantially less energy deposition, temperature differentials of 10°, suggesting that MRTI guidance under laser therapy could provide superior control over the bioheat transfer.
The goal of this work is to deliver a computational model of bioheat transfer that employs real-time, patient specific data and provides real-time high fidelity predictions to be used concomitantly by the surgeon in the laser treatment process. The model employs an adaptive hp-finite element approximation of the nonlinear parabolic Pennes equation and uses adjoint-based algorithms for inverse analysis, model calibration, and adaptive control of cell damage. The target disease of this research are localized adenocarcinomas of the breast, prostate, cerebrum, and other tissues in which a well-defined tumor may form. The algorithms developed also provide a potentially viable option to treat other parts of the anatomy in patients with more advanced and aggressive forms of cancer who have reached their limit of radiation and chemotherapy treatment.
The adaptive data driven computational control system for controlling laser treatment of prostate cancer described is an example of a Dynamic-Data-Driven Applications System (DDDAS) in which simulation models interact with measurement devices and assimilate data over a computational grid for the purpose of producing high-fidelity predictions of physical events.
The development stages uses canines as the host of prostate and brain tumors. The actual laboratory at M.D. Anderson Cancer Center in Houston, TX is connected through a computational grid to the computing center in Austin, TX. Prior to treatment, MRI data is used to generate a finite element mesh of the patient-specific biological domain. The mesh is then optimized to a particular quantity of interest using goal-oriented estimation and adaption. The tool then proceeds to solve an optimal control problem, wherein the laser parameters (location of optical fiber, laser power, etc.) are controlled to eliminate/sensitize cancer cells, minimize damage to healthy cells, and control Heat Shock Protein (HSP) expression. Upon initiation of the treatment process, real-time MRI data is employed to co-register the computational domain and MRTI data is used to calibrate the bioheat transfer model to the biological tissue values of the patient. As new data is given intermittently, computation is compared to the measurements of the real-time treatment and an appropriate course of action is chosen according to the differences seen. A parallel computing paradigm is used to meet the demands of rapid calibration and adapting the computational mesh and models to control approximation and modeling error. From a computational point of view, the orchestration of a successful laser treatment is to solve the problems of co-registration, calibration, optimal control, and mesh refinement invisibly to the surgeon, and merely provide the surgeon with an interface to the optimal laser parameters during treatment. A schematic of the control loop is shown in Figure 1. The mathematical characterization of so-called HSP expression is also described as are laboratory procedures for measuring tissue damage, mesh generation, and the development of the computational infrastructure for calibration procedures, verification and validation procedures, and inverse modeling and sensitivity analysis.
The first step in simulation is to construct a finite element mesh used for the governing bioheat transfer equations. The mesh generation involves two steps: (1) construction, from MRI data, of a finite element mesh to represent the geometry and (2) overlaying the MRTI temperature field onto the finite element mesh.
As an initial test of the imaging system when employed in experiments with laboratory animals, we first consider data pertaining to a mouse. A sample of the MRI data used to construct a mesh of the mouse and tumor is shown in Figure 3. In the case of prostate cancer, prostate tumor cells were inoculated in the hind legs of a mouse and grown to a tumor burden of less than 1.0 cc, Figure 2. The hexahedral mesh of the tumor (blue) and tissue (yellow), shown in Figure 2, was created by a semi-automatic segmentation method adapted to find the interface boundaries of tumor and tissue. Cubic spline and lofting methods are applied to obtain smooth boundaries from the segmented MRI data [16, 21].
Spatio-temporal temperature distribution is measured during the laser treatment with update times less than 6 seconds per image and thickness between planes of 3.5 mm. The temperature field, measured in the biological domain, is first filtered to eliminate any noise and then nodal temperature values are assigned to the mesh by taking the interpolant of the MRTI temperature data.
Heat shock proteins are a latent defense mechanism built into all cells, including cancer cells, that provide enhanced viability of tumor tissue indirectly by preventing a damaged cell from going into apoptosis, resulting in the recurrence of cancer. Knowledge of temperature history versus time during treatment has been used to predict thermal necrosis in regions where damage is severe, but in regions where temperatures are insufficient to coagulate proteins, the results and subsequent effects have been difficult to predict. This is due, in part, to the expression of heat shock protein (HSP) in the regions of thermal stress. Consequently, knowledge of the thermal dose necessary to activate or de-activate HSP expression as a function of temperature and time in the affected tissue [16, 15] can be critical in planning and implementing an effective thermal treatment by laser surgery.
Heat Shock Proteins are a family of gene products expressed in higher concentrations in the presence of environmental stresses. The name vastly understates the HSP family’s astounding versatility. These proteins have been identified as critical components of cell survival under adverse environmental conditions. Various levels of HSP species indicate the health or likelihood of cell proliferation or drug resistance. Also, measures of cell damage as a function of temperature and time for a specific patient is a critical indication of the effectiveness of thermo-therapies. Knowledge of the temperatures necessary to elicit interaction in tumor resistance and cell damage is essential to effectively produce a desired tissue response and surgical outcome. The work of Rylander  developed a model for HSP expression and cell damage based on an Arrhenius model for a mouse. Figure 4 shows the comparison of the experimentally determined HSP expression with the predicted values of (1). An empirical formula for the concentration of HSP expression, , at temperature, u, and time, t, based on laboratory tests is :
where α, β, and γ are time independent parameters that may depend on temperature, with γ > 1.
Cellular damage is measured in terms of the damage fraction FD predicted by means of an Arrhenius integral formulation .
where C0 is the initial concentration of healthy cells, Ct the concentration of healthy cells after heating at time t, A the pre-exponential scaling factor, Ea the activation energy of the injury process, R the universal gas constant, and u the absolute temperature.
Driving the prediction of the HSP and cellular damage is the temperature field produced by the well-known Pennes bioheat transfer model . Liu  has shown Pennes model  to give good results for prediction of temperature field in the prostate. The formal statement of Pennes bioheat model is as follows
Find the spatially and temporally varying temperature field u (x, t) such that
given the Cauchy and Neumann boundary conditions
and the initial condition
Here and are bounded functions of u, determined empirically to be of the form
cp and cblood are the specific heats, ua the arterial temperature, ρ is the density, and h is the coefficient of cooling. The laser source term Qlaser is a linear function of laser power and exhibits exponential decay with distance from the source.
P = is the laser power, μa, μs = are laser coefficients related to laser wavelength and give probability of absorption of photons by tissue, γ is the anisotropy factor, and x0 = is the position of laser photon source. Tables 1 and and22 contain constitutive data from the work of Rylander .
Our weak form of the Pennes bioheat transfer model is as follows:
Given a set of model, β, and laser, η, parameters,
Find u (x, t) H1([0, T], H1(Ω)) s.t.
B(u, β; v) = F(η; v) ∀v
where the explicit functional dependence on the model parameters, β = (k0, k1, 3, 3, ω0(x), ω1, 3, 3), and laser parameters, η= (P(t), x0, μa, μs), expressed as follows
The control loop involves three major problems: calibration of the Pennes model to patient specific MRTI data, optimal positioning and power supply of the laser heat source, and computing goal oriented error estimates. Each of these problems is formally stated below.
The problem of model calibration is to find the set of thermal conductivity parameters k0, k1, 3, 3, and blood perfusivity parameters ω0(x), ω1, 3, 3 that minimize the norm of the difference between the predicted temperature field and the experimentally determined temperature field at each time instance of the experimental data. ω0(x) is allowed to vary over the spatial dimension as the blood perfusivity within the necrotic core of a cancerous tumor is expected to be significantly lower than the surrounding healthy tissue. The calibration problem may be stated as follows:
Given a set of laser parameters, η 0, and an experimentally determined temperature field at time instances tn, n = 1, 2, … Nexp within the region Ωχ Ω
find the best combination of model coefficients, β* , that produces the temperature field, u* , such that
The problem of optimal laser heating control is to find the position, x0, power, P (t), and laser frequency μa, μs that maximizes damage to the cancerous tissue while minimizing damage to healthy tissue in some metric. The metric shown here is the 2 ([0, T ]; 2(Ω)) norm of the difference between the computed solution and an ’ideal’ temperature field. Other metrics may be the 2(Ω) norm of the difference between the computed damage fraction, equation (2), and an ideal damage field within the biological domain at the end of the treatment process or the 2 ([0, T]; 2(Ω)) norm of the difference between the computed HSP expression field, equation (1), and an ideal HSP expression field.
Optimal control of the laser source term is stated as:
Given a set of model coefficients, β 0, and the ideal temperature field that maximizes damage to cancerous tissue while minimizing damage to healthy tissue
where ΩH and ΩC are the domains of the healthy and cancerous tissue respectively, find the best combination of laser position and power, η* , that produces the temperature field, u* , such that
where is the appropriate space for the bioheat transfer model.
Goal oriented error estimates  are calculated under the following framework:
We presume that the solution uhp is a more accurate representation of the exact solution u than uh ; likewise an arbitrary quantity of interest, Q, evaluated at Q(uhp) is a more accurate representation of Q(u) than the quantity of interest evaluated at Q(uh). As an explicit example, consider the quantity of interest to be an averaged value of the temperature in a domain defined through the characteristic function, χ(x),
Using the following Taylor series approximations,
we may solve for the adjoint variable php such that
Where (t) is defined as
This quantity may be used to determine when greater time stepping accuracy is needed. (x), (x), and (x) defined as
may be used to determine when spatial mesh refinement is needed.
The calibration and optimal control problems in the control loop may be posed as the following optimization problem
Find q* s.t.
where q is a parameter in the control space ,
and u is the state variable determined by a variational PDE of the form
in the appropriate Hilbert space, .
The demands of real-time computation will only permit the computation of the following Gateaux derivative, δQ, of the objective function.
The gradient is computed by the adjoint method under the assumptions that the solution, objective function, and constraints admit the following Taylor series expansions.
where the Gateaux derivatives are defined below.
Given q , (and hence u by definition), the derivative of the objective function is
δQ(q,) = δqQ(u(q),q; ) − δqC(u(q),q; , p)
where p is the solution to the adjoint problem.
Find p :
δuC(u(q),q; û;, p) = δuQ(u(q),q; û) ∀û
The solution to the adjoint problem, p, implies that
from the first variation of the PDE constraints we have that
the solution follows
The following Taylor series expansion are used in computing the Gateaux derivatives of the variational form of Pennes model:
Notice that the variational form of the adjoint problem for the calibration, optimal control, and goal oriented error estimates are all the same:
where the only difference is source term δuQ(u(q), q; û) which depends on the quantity of interest. The strong form of the adjoint formulation is as follows:
Given the spatially and temporally varying temperature field u (x, t) find the Lagrange multiplier p(x, t) such that
given the boundary conditions
and the terminal condition
When approximations of the temperature field, u(x, t), and the adjoint variable, p(x, t), are known, the first variation of the quantity of interest for the calibration problem may be obtained explicitly as follows.
The first variation for the quantity of interest for optimal control is as follows:
Two main results are presented in this section. First, the overall computational feasibility of developing a laser treatment paradigm in which high performance computers control the bioheat data transferred from a remote site is demonstrated. Second, reliable computational predictions are shown when comparing results obtained using the Pennes model to the experimental MRTI data in canine cerebral tissue.
The typical time duration of a laser treatment is about five minutes. During a five minute span, one set of MRTI data is acquired every 6 seconds. The size of each set of MRTI data is ≈330kB (256×256×5 voxels) and the bandwidth of a commercial T1 Internet connection is ≈300kB/s. Accounting for connection latency, each set of MRTI data can be transferred between Houston and Austin in under two seconds. A finite element mesh consisting of approximately 10000 degrees of freedom is sufficient to resolve the geometry of the biological domain. The execution time of a representative 10 second bioheat transfer simulation is presented in Figure 6. The execution times represent 10 nonlinear state solves of Pennes model (10 one second time steps) combined with 10 linear adjoint solves to be used for calibration, optimal control, and/or error estimates. Computations were done at the Texas Advanced Computing Center on Dual-Core Linux Cluster. Each node of the cluster contains two Xeon Intel Duo-Core 64-bit processors (4 cores in all) on a single board, as an SMP unit. The core frequency is 2.66GHz and supports 4 floating-point operations per clock period. Each node contains 8GB of memory. The fastest time recorded is .67 seconds, meaning that in a real time 10 second span Pennes model can predict out to more than two minutes! Equivalently, in a 10 second times span roughly 14 corrections can be made to calibrate the model coefficients or optimize the laser parameters.
Preliminary computations comparing the predictions of Pennes model to experimental MRTI taken from a canine brain show very good agreement. A manual craniotomy of a canine skull was preformed to allow insertion of an interstitial laser fiber. Thirty-six two dimensional 256×256 pixels MRI images, Figure 7, of the canine brain were acquired. The field view was 200mm × 200mm with each image spaced 1mm apart. A finite element mesh of the biological domain generated from the MRI data is shown in Figure 8. The mesh consists of 8820 linear elements with a total of 9872 degrees of freedom. MRTI thermal imaging data was acquired in the form of five two dimensional 256×256 pixel images every six seconds for 120 time steps. The spacing between images was 3.5mm. The MRTI data was filtered then projected onto the finite element mesh. Figure 9 shows a cutline comparison between the MRTI data and the predictions of Pennes model. It is observed that the results delivered by the computational Pennes model slightly over diffuses the heat profile peaks compared to measured values. However, at early times the maximum temperature value is within 5% of the MRTI value.
Results indicate that adaptive hp-finite element technology implemented on parallel computer architectures and employing modern data transfer infrastructure and thermal imaging modalities are capable of predicting and guiding computer controlled temperature field within a biological domain at very good accuracies. Reliable finite element model simulations of laser therapies can be computed in a fraction of the time that the actual therapy takes place. These prediction capabilities combined with an understanding of HSP kinetics and damage mechanisms at the cellular and tissue levels due to thermal stress, together with recent advances for controlling discretization and modeling errors through adaptive control strategies combined with the modern methods of inverse analysis, calibration, uncertainty quantification, sensitivity analysis, and optimization [22, 12, 11, 14, 2, 8, 10, 9, 1, 3] provide a powerful methodology for planning and optimizing the delivery of thermo-therapy for cancer treatments. Moreover, from a computational point of view, changing the thermal delivery modality merely amounts to changing the source term in the governing PDE. This technology has the potential to be extended to many areas of thermal treatment including RF, microwave, ultrasound, and even cryotherapy applicators.
The support of this project by the National Science Foundation under grant CNS-0540033 is gratefully acknowledged.