|Home | About | Journals | Submit | Contact Us | Français|
Predicting the outcome of thermotherapies in cancer treatment requires an accurate characterization of the bioheat transfer processes in soft tissues. Due to the biological and structural complexity of tumor (soft tissue) composition and vasculature, it is often very difficult to obtain reliable tissue properties that is one of the key factors for the accurate treatment outcome prediction. Efficient algorithms employing in vivo thermal measurements to determine heterogeneous thermal tissues properties in conjunction with a detailed sensitivity analysis can produce essential information for model development and optimal control. The goals of this paper are to present a general formulation of the bioheat transfer equation for heterogeneous soft tissues, review models and algorithms developed for cell damage, heat shock proteins, and soft tissues with nanoparticle inclusion, and demonstrate an overall computational strategy for developing a laser treatment framework with the ability to perform real-time robust calibrations and optimal control. This computational strategy can be applied to other thermotherapies using the heat source such as radio frequency or high intensity focused ultrasound.
Thermotherapy using heat sources such as laser, radio frequency, focused ultrasound, and microwaves provide an alternative to surgery or chemotherapy in the treatment of solid tumors embedded in vital regions. It offers many advantages for being minimally invasive and allowing for targeted treatment [20,22,27,41,42]. Technology advances, such as actively cooled applicators, multi-functional nanoparticles and image-guided probe placement, have made thermotherapy a viable and attractive therapeutic modality [12,16,32,33]. One of the requirements for effective treatment of small, poorly-defined metastases or other tumors embedded in vital organs is to control so-called specificity, which could potentially be achieved, for example, through effective integration of nanoparticles into the laser therapy. In general a successful thermotherapeutic operation requires pre-treatment planning, image-guided surgical control, and post-treatment monitoring. Of key importance to such a procedure is the ability to predict accurately the temperature field and identify appropriately the damaged tissue regions. One of the major challenges associated with predicting the temperature field is that the tissue properties can be difficult to measure and may vary over time during the treatment due to biological alteration in the tissues. Moreover, the complexity of tissue composition and structures makes it difficult to predict the treatment outcome accurately without a reliable mathematical and computational model to characterize the bioheat transfer process.
In this paper, we first give a brief review of bioheat transfer models in the literature with the major parameters being identified. This is followed by a systematic sensitivity study of temperature and corresponding tissue damage with respect to a selected parameter set for both linear and nonlinear bioheat transfer models. In Section 3, a general nonlinear heterogeneous bioheat transfer model is formulated, which is capable of capturing heterogeneity of thermal properties due to anatomic structure difference (i.e. smooth tissue, veins, and arteries), injection of nanoparticles, and the variation in time due to temperature dependent perfusion, provided in vivo temperature data are available. This nonlinear heterogeneous model can be implemented efficiently on a high performance computer system using parallel algorithms, which permits real-time monitoring and control of a thermal surgical process . Numerical examples are presented using a three-dimensional finite element method for both phantom and canine prostate models with an isothermal laser heating source. Results show that this model can reliably characterize changes in tissue properties and accurately predict temperature fields comparable to that measured with an in vivo magnetic resonance temperature imaging (MRTI) technique.
The bioheat transfer processes in living tissue involve conduction, convection, perfusion, metabolism, evaporation, and radiation. The influence of the circulating blood flow through embedded vascular structures arguably plays a key role in the case of soft tissues perfused with blood. In a human body, the heart generates the pressure driving the blood through the branching system of vessels from aorta, arteries, to arterioles, until they reach the level of capillaries. Then, the blood drains into the small venous vessels, flows through larger veins and venae cavae, and eventually returns back to the heart. In this circulation process, there are not only differences in size and structure of these vessels but also in blood temperature. Further distinctions can be made between healthy and cancerous tissues in terms of morphology, vascular structure, and temperature.
To model bioheat transfer processes, the initial work by Pennes  published in 1948 provides a basis for describing diffusion and perfusion with the assumption that no heat transfer between the blood vessels and tissue takes place until the level of capillaries is reached. The heat exchange, in this model, is treated as immediate and complete at this level. Other key assumptions are that the tissue is a homogeneous media, and the thermal energy transfers through perfusion, acting as volumetric heating source or sink proportionally to the difference between local tissue temperature and the blood temperature (considered to be constant). The model neither allows for mass transport through blood vessels nor accounts for nonlinear effects due to the vascular structure in tissues.
Over the last six decades, the Pennes model has been revisited and improved by many investigators (e.g., [5,6,39]). In , Chen and Holmes showed that the assumption regarding heat exchange occurring only at the capillary level is flawed since thermal equilibrium is reached in arterioles with diameter ~60 μm. To improve Pennes model, these investigators considered tissues consisting of solid and fluid parts with the energy being conserved in each part but not allowing for mass transfer between them. Since this model does not take into account the vascular structures in the tissues, the effect of thermal energy transferred via a counter-current produced between veins and arteries is missing. To that end, many attempts have been made (e.g., [1,3,4,17,18,28,29,35,36,38,40,43,44]) by taking into account details of the vascular structure and its nonlinear effect. As a result, bioheat transfer models have become more complex. Moreover, detailed accounts of the vasculature, in addition to unmeasurable material coefficients, need to be addressed to implement these models. For a comprehensive review of bioheat transfer models, we refer to Diller et al. .
Although the Pennes model has been shown to be limited due to its overly simplified perfusion term , it remains a usable model merited by its simplicity and ease of implementation. Moreover, the blood perfusion term in the Pennes equation serves as an adjustable parameter that can be chosen to provide reasonable agreement between the model predictions and experimental measurements .
In this section, we present a general form of the bioheat equation for heterogeneous soft tissues [24,25]. Let Ω be a bounded domain in 3 with a Lipschitz continuous boundary, i.e. viable tissue with a tumor. Eq. (1) is based on the Pennes bioheat equation . However, we permit the conductivity k and perfusion coefficient ω to be spatially heterogeneous and temperature dependent
The specific heat of the tissue and blood are given by cp and cb, respectively, ua is the arterial temperature, ρ is the density of the tissue, and Qlaser(β,x, t) represents a heating source per unit volume. Based on observations of the thermal response, the thermal conductivity, k(u,x,β), and blood perfusion, ω(u,x, β), which are assumed to be bounded functions of the temperature field, u = u(x, t), where
and k0(x), k1, k2, k3, ω0(x), ω1, ω2, and ω3 are parameters of the diffusivity and perfusion coefficient functions defined above. The heterogeneity of the tissue is modeled in the first terms in the thermal conductivity and perfusion. The source term for laser heating is assumed to take the form
where P is the power of the laser, x0 is the position of the laser tip, μa and μs, respectively, are the absorption and scattering coefficients, g is the anisotropic factor, μtr = μa + (1 − g)μs and (see Section 5.1 for discussion of these parameters).
On the Cauchy boundary,
where the convection coefficient is denoted as h and u∞ represents the ambient temperature. The parameter is the prescribed heat flux on the Neumann boundary
The temperature field is propagated forward in time from a given initial condition, u0
Next, we discuss the optimization algorithms with all model parameters for the perfusion, thermal conductivity, and laser source term are included in an array β defined as
For predicted treatment outcome to be clinically relevant, it is necessary to formulate the mathematical setting of the thermotherapeutic protocol based on cell and tissue damages. It is too often the case that the temperature is used to formulate the objective functional, although the temperature is the control variable that can be measured in vivo. However, it cannot be overemphasized that cell and tissue damages are the ultimate measure whether a thermotherapy is successful. An optimized temperature field does not necessarily imply an optimized tissue damage [31,32]. Therefore, it is very important to formulate a biologically relevant objective functional.
The thermal injury resulting from a heat source can be defined by the damage field, D(x), which is also dependent on the temperature history. Let ϕ denote the functional that acts on the temperature field and returns the corresponding damage field
where u denotes temperature and τ the final heating time. Assume there exists an ideal damage field, Dideal(x), that maximizes damage to the cancerous region while maintaining functionality of the healthy region. As a concrete example, consider the ideal damage field to be a discontinuous function of x Ω. Within the cancerous region there is a threshold damage value, denoted as , above which all cancerous cells will be killed. Similarly, within the healthy region there is a threshold damage value, denoted as , below which the functionality is maintained for all healthy cells
To determine the effectiveness of a laser treatment, we define the damage field resulting from the application of a laser heat source ϕ(u(β,x, t)) = D(x) satisfies the ideal conditions:
where |·| denotes the Lebesgue measure of the set, and
As mentioned previously, the goal of the treatment is to maximize the damaged volume of cancerous tissues that have reached a damage threshold to eradicate the cancerous cells and minimize the volume of healthy tissue that reach a threshold beyond which the healthy tissue begins to lose functionality.
Since the temperature field is a continuous function, the resulting damage field is expected to be a continuous function. Therefore, Eq. (2) is unattainable. There will always exist a transitional region where the damage field does not satisfy the ideal criteria. In this region, healthy cells will lose functionality and/or cancerous cells will not be destroyed. An one-dimensional illustration is shown in Fig. 1. These transitional regions are highlighted within shaded blue boxes.
Therefore, we seek the optimal treatment that can be achieved in order to minimize the region in which the ideal criteria is not achieved. In other words, we look for solutions to minimize the L2-norm of the difference, denoted as Q, between the ideal damage field and the computed damage field
Thus, an optimization problem can be formulated to minimize above objective functional over the parameter space as follows:
However, it is possible that a completely acceptable solution to Eq. (2) may have a large L2 error (ref. Fig. 1). Therefore, during the solution process of Eq. (3), the ratio of the volume of the domain that satisfies the ideal criteria to the total volume of the domain should be used as a measure of the degree of laser treatment effectiveness to which Eq. (3) satisfies Eq. (2),
In the event that the ratios (4) resulting from the solution to (3) are unsatisfactory, the criteria may be relaxed by allowing certain percentage of the healthy cells to be destroyed along with the cancerous cells, i.e.
where ε is a small region around the boundary of ΩC and ε ΩH. In clinical situations, the parameter ε is the margin area determined by experienced surgeons.
In a previous study , we developed a two-state model to quantify the damage field from the thermal injury in order to overcome the inaccuracy and numerical instability of the classical Arrhenius model, which is an empirical one rooted from characterization of kinetics for chemical reactions. It is not until late 1940s that this empirical law was adopted to describe tissue damage [14,15,21]. The two-state model , on the other hand, was developed based on statistical thermodynamic principles. For convenience, let ϕ represent the damage field for both the Arrhenius and the two-state models
In the Arrhenius model, the parameters Ea, A, and R are known constants that represent the activation energy, the frequency factor, and the universal gas constant respectively. In the two-state model, is an activation energy function that depends on the temperature u and time t, where h, α, and β are known constants determined by in vitro cellular experiments. The parameter κ is the Boltzmann constant. The upper limit τ in the integrals is the final time of heating. It is clear that the numerical values of and for the Arrhenius model will be different from and for the two-state model.
Other relevant biomarkers (quantities of interest) related to tumor recurrence are heat shock proteins (HSPs). After the thermotherapy, HSP levels usually increase in both healthy and cancerous cells. Previous studies have documented the link between HSP27 and HSP70 (number denotes molecular weight in kilodaltons) and tumor recurrence due to their roles in multi-drug resistance and regulation of apoptosis through modulation of p53 oncogene. In human prostate, breast, and cervical cancers, over-expression of HSP27 and HSP70 are poor prognostic biomarkers (e.g., [7,19]). Increased HSP27 and HSP70 levels during thermal treatment, by exposure to sub-lethal temperatures, are very likely to result in surviving of tumor cells via inhibition of apoptosis and increased resistance to subsequent radio frequency and chemotherapies. Prediction of the thermally induced HSP27 and HSP70 levels following the thermotherapy can provide an indication of the likelihood of tumor recurrence based on the measured thermal dose, which can be included in the objective functional [31,33]. Tumor regions having HSP expression above the basal level (1 mg/ml) are insufficiently treated with the thermal dose. It could lead to a high likelihood of tumor re-growth and resistance to subsequent therapies thereby increasing the likelihood of tumor recurrence. Tumor regions with extensive temperature elevation, high levels of injury, and minimal HSP expression (less than 1 mg/ml) due to denaturation of proteins will confirm the effectiveness of a lethal thermal dose to the targeted tissue regions. In a previous experimental study , we proposed that the HSP level for both healthy and cancerous prostate cells (PC3 and RWPE-1), denoted as H(u, t), can be formulated using the following relationship:
where H represents the level of HSP expressed as a function in terms of time t and temperature u. Parameters α, β, and γ (>1) are to be determined by cellular experiments. It is believed that this is the first quantitative HSP expression model with a given thermal dose based on measured thermally induced HSP expression in prostate cells. This model was incorporated into the treatment planning and optimization model for control of the clinical outcome [31–33].
Now, we define the mathematical problem with partial differential equation (PDE) constrains based on the general form of the nonlinear bioheat transfer Eq. (1).
Let L2(0, τ;H1(Ω)) be the trial space and the dual space of , where τ the final time. The problem can be formally stated as follows:
Given a set of model parameters, β, find
where : → . Denoting the duality pairing by ·,· × , we have
Based on the theoretical results from Showalter  on the space L2(0,τ;H1(Ω)), along with time derivatives in (L2(0,τ;H1(Ω)))′, there exists a pointwise solution with respect to time as long as the initial condition u(x,0) = u0 is well-defined.
Next, we determine the parameter space, over which the optimization is performed, by defining
To incorporate the damage model-based optimization with the PDE constraint, the variational framework must be posed in the context of a minimization problem for the objective functional, Q, over parameter space,
Find β* s.t.
where u is the state variable determined by a variational form of the PDE,
that defines the governing PDE, which is nonlinear with respect to the first two arguments, u, β, but linear with respect to the test function, ν. In the case of the bioheat transfer equation, Eq. (1), it take the following form:
where u, ν L2([0,τ],H1(Ω)) and the control variable, β, belongs to the relevant parameter space, β , as in Eq. (9).
Due to the stringent time constraints for real-time computation, we choose not to use the Hessian of the objective function to compute the rate of change with respect to each model parameter. Instead, we compute the gradient of the objective function using quasi-Newton method  to solve the optimization problem, in which only limited-memory variable metric is needed. A detailed account of the algorithms for the gradient computation can be found in [8,25].
It is important to note that the efficiency and accuracy of model calibration are very important for reliable real-time surgical control. In order to achieve high accuracy, it is necessary to perform model calibration to capture patient specific tissue properties as well as dynamic changes due to coagulation and/or nanoparticles inclusion, which is often used to enhance the effectiveness of thermotherapies [10,12,16]. In other words, these factors are the sources of tissue heterogeneity. The general idea is to use in vivo MRTI data and efficient inverse analysis algorithms to solve for tissue properties that is spatial dependent and change over the time.
The primary quantities of interest in our study for characterizing soft tissues are the thermal conductivity k and the perfusion coefficient ω. Although the thermal conductivity may vary little with respect to temperature, it is believed that the perfusion coefficient is highly temperature and time-dependent. Using optimization framework, we are able to solve for both thermal conductivity and perfusion coefficients in real-time, which captures their spatial and temporal dependence. The main idea is to solve for tissue properties using inverse analysis formulated by a PDE-constrained optimization problem, provided that measured temperature filed (MRTI data) are given. We refer this process as calibration.
The objective functional for the calibration problem is defined as the difference, in space–time norm, between the computed temperature field, u(x, t), and the temperature field measured by the MRTI experiments, uexp(x, t) integrated over the entire biological domain Ω and time duration of interest [0,τ].
The demands of real-time prediction place stringent demands on the speed and accuracy of the algorithms used to calibrate the model . Efficient algorithms to calibrate the generalized nonlinear bioheat transfer model have been studied in . In this section, we address the following important questions: Will the optimizer converge to a solution in the time allowed for real-time computation? Furthermore, is the convergence of the optimizer robust? “Robust” in this context implies that the optimization process will converge independently of the finite element discretization, independent of the initial guess of the model parameters in physical parameter space, and independent of which model parameters are optimized. With regard to the latter criteria, it is observed that the absorption coefficient μa, the first term in the thermal conductivity, k0, and the first term in the blood perfusion, ω0 are the most sensitive parameters in the calibration problem. However, it is also found that the calibration process with only two parameters (k0 and ω0) is insuf.- cient to minimize the objective functional. It seems that three parameters (μa, k0, and ω0) are the minimum set for the β array for the calibration problem.
On the other hand, the experimental system for model calibration was based on phantom experiments. The phantom samples were made of an excised canine prostate embedded within a 1% agar gel. Using interstitial laser heating with real-time temperature feedback measured by MRTI, it was demonstrated that the calibration algorithm was robust and efficient for real-time control. The user interface, shown in Fig. 2, illustrates the phantom geometry, predicted temperature field, predicted damage field, cutlines of the thermal images, cutlines of the heating predicted by Pennes model, and the planned laser power as a function of time.
Calibration was conducted using a time series of thermal images taken throughout the heating of the phantom material. Permutations over the highest and lowest values of physically meaningful model parameters as the initial guess for calibration performed. Finite element meshes of 33,096 degrees of freedom and 87,830 degrees of freedom were used in the study. It is obvious that the convergence properties depends on the choice of the initial model parameters It was found that the optimization problem converged under 10 function/gradient evaluations for most of the scenarios considered. Previous results  indicate that it requires approximately 10 function/gradient evaluations for a real-time computation during a typical laser surgery for cancer treatment of mid-size tumors. To achieve the fastest convergence, in our case, the initial model parameters should be chosen to be close to that in the linear model. In that case, the optimization algorithm will solve for the nonlinearity efficiently.
The introduction of nanoparticles into tumor regions can greatly enhance the effectiveness of thermotherapy. However, these nanoparticles inside the tissue can alter the original tissue properties. For example, nanoshell-mediated laser surgery is able to regulate thermal energy into target regions delivered by optical fibers to provide a lethal dose of heat while minimizing damage to surrounding tissue . A nanoshell is composed of a spherical core of a particular compound surrounded by a shell of a few nanometer in thickness. In our study, the nanoshell has a gold shell with a silicon dioxide core, which was manufactured by Nanospectra Biosciences, Inc. Nanoshells can act as intense infrared absorbers which increase the thermal deposition of laser energy into tumors. In particular, nanoshells provide a potential means to enhance the delivery of laser-induced thermal energy via distributing the heat source from the fiber to the surrounding vasculature and/or provide a highly conformal and targeted approach to laser-induced thermal therapy in which normal tissue is spared and tumor tissue is ablated with a high level of specificity. Previous studies [12,13] show that the inclusion of nanoparticles can be used as a pragmatic strategy for cancer treatment. The algorithm discussed previously can serve the same purpose to estimate the effective tissue properties due to nanoparticle inclusion , by performing inverse analysis using the general optimization framework discussed in Section 3.4.
A crucial element in a high fidelity bioheat transfer model is capturing the heating source accurately. In the case of laser-induced thermotherapy, this requires the laser fluence in the tissue to be modeled correctly. An analytic solution of the transport equation, solved with a spherically isotropic point source, has been used for the fluence in this study.
where P represents the power of laser in Watts as a function of time, μtr = μa + (1 − g)μs, and (see Chapter 6 in  for details).
Three major parameters in this model are μa, μs, and g, which represent the average number of photons that are absorbed or scattered per unit length, and the expected value of the cosine of the scattering angle, respectively. Intuitively, values of g close to 1 indicate a medium which scatters light in a more of a forward direction while values close to −1 indicate a medium which generally scatters light backwards. Soft tissue is generally a forward scattering medium with values of g greater than 0. Each parameter, when considered for living tissue, should be a function of space, light wavelength, and temperature. As a first order approximation, these parameters are often taken to be functions of wavelength only. Under this approximation, μa, μs, and g are each constants for the current application, since light emitted from a laser is of a single wavelength.
Experimental results indicate that changes in temperature may result in significant changes in these optical parameter values . This raises the question of whether such a change in the heating source parameters affects the overall thermal profile enough to warrant modeling them as functions of temperature. To that end, a sensitivity study was conducted, the results of which were published in  and are summarized here.
Various simulations were performed for each heating protocol of a tumor on the back of a mouse. Initial values of parameters were taken as μa = 3.14 (1/cm), μs = 14.2 (1/cm) and g = 0.7.
Fig. 3 shows representative snap shots of temperature profiles for different values of μa with initial values for μs and g each after 3 min and 20 s of heating, a typical duration time of thermal therapies.
Generally, higher values of μa and μs were associated with an increase in heating, whereas an increase in g was associated with a decrease in heating. The thermal profile appeared less sensitive relative to μa and μs as opposed to g. The results indicated that characterizing μa and μs with constants is justifiable for this form of the fluence, however, further investigation should be done to determine the amount of change g is expected to undergo during the treatment.
The Arrhenius damage model-based study was performed for the scenario depicted in Figs. 4 and and55 based on a finite element model for a canine prostate. A series of four-minute treatments were simulated. The power was varied from 3.0–6.0W in increments of 0.5 W. Values at the physical upper and lower bounds over the range of the thermal conductivity, perfusion, absorption coefficient, and laser power were considered. Arrhenius damage model data for in vitro cellular were adopted from Rylander . As expected, the large coefficients of the Arrhenius model, A ≈ 1017, are observed to make the computation of the damage model-based objective functional very sensitive to the model parameters. A small change in the temperature field produces a large change in the damage. Furthermore, the gradient is proportional to the large coefficients of the Arrhenius model. The optimization process is expected to change the initial guess in parameters space proportional to the gradient. The result being that the function evaluation is expected to require a significant number of evaluations to converge to a solution and ultimately unable to make a real-time Arrhenius damage model optimization.
The mathematical models and computational methods described in this paper are at the core of a new computational paradigm for predicting the outcome of patient specific, thermal therapy in prostate cancer treatment. Key to this effort are high fidelity predictions of the temperature field, cell damage, and HSP expression, all of which mandate the use of robust and accurate models along with quality data for the tissue perfusion, the laser source term, and heterogeneous, possibly time-dependent, soft tissue properties. In this regard, the work presented herein has taken preliminary steps toward investigating the use of various classes of models and has performed a sequence of parameter studies to determine the relative sensitivity of these models to the major parameters.
Based on the sensitivity tests reported in this study, the following key results are noted:
In conclusion, we have shown that high performance computing coupling with real-time in vivo thermometric measurement using magnetic resonance temperature imaging (MRTI), we are able to calibrate bioheat transfer models in terms of both nonlinearity and heterogeneity resulting from vasculature, tissue damage, or nanoparticle inclusions using a unified computational framework in order to make reliable predictions for meeting the goal of patient-specific treatment planning and real-time thermotherapeutic control.
As Prof. Oden’s current and former students, and former post-doctoral fellow, we are deeply grateful for his encouragement, guidance, leadership and active involvement in this new endeavor of computational cancer research. Also, we would like to acknowledge that the work reviewed in this article reflects a fruitful collaboration of a multidisciplinary team, which is led by Prof. Oden and consists of K.R. Diller, J. Hazle, J.C. Browne, I. Babuška, L. Bidaut, L. Demkowicz, C. Bajaj, A. Elliott, J. Zhang, S. Goswami, S. Khoshnevis, S. Prudhomme, R.J. Stafford, and A. Shetty, from the University of Texas at Austin and M.D. Anderson Cancer Center in Houston. Finally, the support of this work by the National Science Foundation under Grant CNS-0540033 and the National Institutes of Health under Grant 7K25CA116291-02 (YF) are gratefully acknowledged.