PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Methods Appl Mech Eng. Author manuscript; available in PMC 2010 May 17.
Published in final edited form as:
Comput Methods Appl Mech Eng. 2009; 198(21): 1742–1750.
doi:  10.1016/j.cma.2008.12.027
PMCID: PMC2871336
NIHMSID: NIHMS196079

Optimization and real-time control for laser treatment of heterogeneous soft tissues

Abstract

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.

Keywords: Bioheat transfer model, Nonlinearity, Parameter estimation, Soft tissue properties, Nanoparticle inclusion, Parallel computing

1. Introduction

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 [25]. 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.

2. Mathematical model of bioheat transfer process

2.1. Review of previous work

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 [26] 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 [6], 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. [9].

Although the Pennes model has been shown to be limited due to its overly simplified perfusion term [6], 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 [35].

2.2. A nonlinear bioheat transfer model

In this section, we present a general form of the bioheat equation for heterogeneous soft tissues [24,25]. Let Ω be a bounded domain in R3 with a Lipschitz continuous boundary, i.e. viable tissue with a tumor. Eq. (1) is based on the Pennes bioheat equation [26]. However, we permit the conductivity k and perfusion coefficient ω to be spatially heterogeneous and temperature dependent

ρcput·(k(u,x,β)u)+ω(u,x,β)cb(uua)=Qlaser(β,x,t)inΩ.
(1)

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

k(u,x,β)=k0(x)+k1atan(k2(uk3)),ω(u,x,β)=ω0(x)+ω1atan(ω2(uω3))

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

Qlaser(β,x,t)=3μaP(t)μtreμeff||xx0||4π||xx0||,

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 μeff=3μaμtr (see Section 5.1 for discussion of these parameters).

On the Cauchy boundary,

k(u,x,β)u·n=h(uu)onΩC,

where the convection coefficient is denoted as h and u represents the ambient temperature. The parameter An external file that holds a picture, illustration, etc.
Object name is nihms196079ig1.jpg is the prescribed heat flux on the Neumann boundary

equation image

The temperature field is propagated forward in time from a given initial condition, u0

u(x,0)=u0inΩ.

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

β(k0(x),k1,k2,k3,ω0(x),ω1,ω2,ω3,P(t),μa,μs,x0).

3. Biologically based optimization

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.

3.1. The 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

ϕ:L2(0,τ,H1(Ω))L2(Ω),ϕ(u(β,x,t))=D(x)xΩ,

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 [set membership] Ω. Within the cancerous region there is a threshold damage value, denoted as DCideal, above which all cancerous cells will be killed. Similarly, within the healthy region there is a threshold damage value, denoted as DHideal, below which the functionality is maintained for all healthy cells

Dideal(x){DCideal,xΩC,DHideal,xΩH.

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:

YC=ΩCandYH=ΩH,
(2)

where |·| denotes the Lebesgue measure of the set, and

YC{xΩC:D(x)DCideal},YH{xΩH:D(x)DHideal}.

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.

Fig. 1
(a) The transitional damage regions (highlighted in shaded blue boxes) that the ideal criteria are not satisfied in the cases of either healthy cells will lose functionality and/or cancerous cells will not be destroyed. (b) a solution that is acceptable ...

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

Q(u(β,x,t))=12||ϕ(u(β,x,t))Dideal(x)||L2(Ω)2=12||D(x)Dideal(x)||L2(Ω)2.

Thus, an optimization problem can be formulated to minimize above objective functional over the parameter space P as follows:

minβPQ(u(β,x,t)).
(3)

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),

YCΩCandYHΩH.
(4)

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.

Dεideal(x){DCideal,xΩC+ε,DHideal,xΩHε,

where ε is a small region around the boundary of ΩC and ε [subset or is implied by] ΩH. In clinical situations, the parameter ε is the margin area determined by experienced surgeons.

3.2. Cell damage models

In a previous study [11], 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 [11], 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

ϕ(u(β,x,t))=D(x)={0τAeEaRudtArrhenius,0τeE2κu1+eE2κudtTwoState.

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, E2=(hu+αt+β) 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 DCideal and DHideal for the Arrhenius model will be different from DCideal and DHideal for the two-state model.

3.3. Characterization of heat shock proteins

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 [32], 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:

H(u,t)t=(αβtγ1)H(u,t),
(5)

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 [3133].

3.4. PDE-constrained optimization

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 An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpg [equivalent] L2(0, τ;H1(Ω)) be the trial space and An external file that holds a picture, illustration, etc.
Object name is nihms196079ig3.jpg the dual space of An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpg, where τ the final time. The problem can be formally stated as follows:

Given a set of model parameters, β, find

equation image
(6)

where An external file that holds a picture, illustration, etc.
Object name is nihms196079ig4.jpg: An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms196079ig3.jpg. Denoting the duality pairing by left angle bracket·,·right angle bracket An external file that holds a picture, illustration, etc.
Object name is nihms196079ig3.jpg × An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpg, we have

equation image
(7)

and

equation image
(8)

Based on the theoretical results from Showalter [34] 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.

It can be readily shown that Eq. (6) has a solution provided that the operator An external file that holds a picture, illustration, etc.
Object name is nihms196079ig4.jpg: An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms196079ig3.jpg is bounded, coercive, and of type x2133 (see, e.g., [34]).

Next, we determine the parameter space, over which the optimization is performed, by defining

P={βL(Ω)×R3×L(Ω)×R3×L([0,τ])×R5:0<k<k0(x)+k1atan(k2(uk3))<k<0<ω<ω0(x)+ω1atan(ω2(uω3))<ω<}.
(9)

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, P

Find β* [set membership] P s.t.

Q(u(β),β)=infβPQ(u(β),β)

where u [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpg is the state variable determined by a variational form of the PDE,

equation image

Using the same notation in Eqs. (7) and (8), we have

equation image

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:

equation image
(10)

where u, ν [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms196079ig2.jpg [equivalent] L2([0,τ],H1(Ω)) and the control variable, β, belongs to the relevant parameter space, β [set membership] P, 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 [2] 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].

4. Real-time model parameter identification

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.

4.1. Estimation of thermal properties of heterogeneous soft tissues

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,τ].

Q(u(β),β)=12||u(β)(x,t)uexp(x,t)||2.

The demands of real-time prediction place stringent demands on the speed and accuracy of the algorithms used to calibrate the model [25]. Efficient algorithms to calibrate the generalized nonlinear bioheat transfer model have been studied in [10]. 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.

Fig. 2
A screen shot of the user interface for the phantom test: the geometry fused with MRTI data (upper-left window), 2D section of predicted temperature field (uppermiddle window, which is also used to show HSP expression), predicted damage field (lower-left ...

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 [25] 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.

4.2. Estimation of effective tissue properties with nanoparticle inclusion

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 [16]. 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 [10], by performing inverse analysis using the general optimization framework discussed in Section 3.4.

5. Parameter sensitivity in bioheat transfer model

5.1. Laser-heating model and parameter sensitivity

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.

Φ(xx0)=3P(t)μtr4π||xx0||eμeff||xx0||,
(11)

where P represents the power of laser in Watts as a function of time, μtr = μa + (1 − g)μs, and μeff=3μaμtr (see Chapter 6 in [37] 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 [23]. 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 [8] 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.

Fig. 3
Images shown here are of the temperature field generated with an isotropic source term after 3 min and 20 s of heating. From left to right, values of μa are taken to be 0.44, 3.14, and 5.0 (1/cm). Values of the other parameters are μs ...

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.

5.2. Damage model sensitivity

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 [30]. 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.

Fig. 4
The Arrhenius damage model objective function is evaluated over parameter space and presented as a semi-log plot. Permutations of the high and low values of physically meaningful model parameters are used. μa, k0, k1, k2, k3, and ω0 are ...
Fig. 5
A screen shot of the canine prostate experiment with damage based optimal control using Arrhenius damage model. The geometry fused with MRTI data (upper-left window), 2D section of predicted temperature field (upper-middle window, which is also used to ...

6. Discussion and conclusions

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:

  • Changes in the absorption and scattering coefficients in the laser fluence do not dramatically affect the thermal profile but the anisotropic factor may. Thus, further investigation is needed to determine the amount of change the anisotropic factor is expected to undergo during laser treatment.
  • The frequency factor A (on the order of 1017) in the Arrhenius damage model makes the computation of the damage model objective function very sensitive to the model parameters. This high sensitivity is expected to result in a significant number of function evaluations to converge to a solution and ultimately obstructs real-time Arrhenius damage model optimization convergence. We will implement two-state model [11], which is less sensitive to model parameter.

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.

Acknowledgments

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.

References

1. Arora D, Skliar M, Cooley D, Blankespoor A, Moellmer J, Roemer R. Nonlinear model predictive thermal dose control of thermal therapies: experimental validation with phantoms. American Control Conference; 2004.
2. Benson SJ, McInnes LC, Moré J, Sarich J. TAO user manual (rev. 1.8), Technical Report ANL/MCS-TM-242. Mathematics and Computer Science Division, Argonne National Laboratory; 2005. < http://www.mcs.anl.gov/tao>.
3. Charny CK, Levin RL. Bioheat transfer in a branching countercurrent network during hyperthermia. J Biomech Engrg. 1989;111:263–270. [PubMed]
4. Chato JC. Heat transfer to blood vessels. ASME J Biomech Engrg. 1980;102:110. [PubMed]
5. Chato JC, Lee RC. The future of biothermal engineering. Ann NY Acad Sci. 1998;858(1):1. [PubMed]
6. Chen MM, Holmes KR. Microvascular contributions in tissue heat transfer. Ann NY Acad Sci. 1980;335:137. [PubMed]
7. Cornford PA, Dodson AR, Parsons KF, Desmond AD, Woolfenden A, Fordham M, Neoptolemos JP, Ke Y, Foster CS. Heat shock protein expression independently predicts clinical outcome in prostate cancer. Cancer Res. 2000;60:7099–7105. [PubMed]
8. Diller KR, Oden JT, Bajaj C, Browne JC, Hazle J, Babuška I, Bass J, Bidaut L, Demkowicz L, Elliott A, Feng Y, Fuentes D, Goswami S, Hawkins A, Khoshnevis S, Kwon B, Prudhomme S, Stafford RJ. Advances in Numerical Heat Transfer. Vol. 3. Taylor & Francis; 2008. Computational infrastructure for the real-time patient-specific treatment of cancer. (Chapter 9)
9. Diller KR, Valvano JW, Pearce JA. Bioheat transfer, The CRC Handbook of Mechanical Engineering. (2) 2005;4:278–357.
10. Feng Y, Fuentes D, Hawkins A, Bass J, Rylander MN, Elliott A, Shetty A, Stafford RJ, Oden JT. Nanoshell-mediated laser surgery simulation for prostate cancer treatment. Special Issue of Comput Bioeng, Engrg Comput. 2008 doi: 10.1007/S00366-008-0109-Y. [PMC free article] [PubMed] [Cross Ref]
11. Feng Y, Oden JT, Rylander MN. A two-state cell damage model under hyperthermic conditions: theory and in vitro experiments. J Biomech Engrg. 2008;130(041016):1–10. [PMC free article] [PubMed]
12. Feng Y, Rylander MN, Bass J, Oden JT, Diller K. Optimal design of laser surgery for cancer treatment through nanoparticle-mediated hyperthermia therapy. NSTI-Nanotech. 2005;1:39–42.
13. Ferrari M. Cancer nanotechnology: opportunities and challenges. Nature Rev Cancer. 2005;5(3):161–171. [PubMed]
14. Henriques FC. Studies of thermal injury V. The predictability and the significance of thermally induced rate processes leading to irreversible epidermal injury. Arch Pathol. 1947;43:489–502. [PubMed]
15. Henriques FC, Moritz AR. Studies of thermal injury. I. The conduction of heat to and through skin and the temperatures attained therein. A theoretical and an experimental investigation. Am J Pathol. 1947;23(4):530–549. [PubMed]
16. Hirsch LR, Stafford RJ, Bankson JA, Sershen SR, Rivera B, Price RE, Hazle JD, Halas NJ, West JL. Nanoshell-mediated near-infrared thermal therapy of tumors under magnetic resonance guidance. Proc Natl Acad Sci. 2003;100(23):13549–13554. [PubMed]
17. Jiji LM, Weinbaum S, Lemons DE. Theory and experiment for the effect of vascular microstructure on surface tissue heat transfer. Part II. Model formulation and solution. ASME J Biomech Engrg. 1984;106:333. [PubMed]
18. Klinger HB. Heat transfer in perfused biological tissue. I. General theory. Bull Math Biol. 1974;36:403. [PubMed]
19. Madersabacher S, Grobl M, Kramer G, Dirnhofer S, Steiner G, Marberger M. Regulation of heat shock protein 27 expression of prostatic cells in response to heat treatment. Prostate. 1998;37:174–181. [PubMed]
20. Masters A, Steger AC, Lees WR, Walmsley KM, Bown SG. Interstitial laser hyperthermia: a new approach for treating liver metastases. Br J Cancer. 1992;66(3):518–522. [PMC free article] [PubMed]
21. Moritz AR, Henriques FC. Studies of thermal injury. II. The relative importance of time and surface temperature in the causation of cutaneous burns. Am J Pathol. 1947;23:695–720. [PubMed]
22. Nagata Y, Hiraoka M, Akuta K, Abe M, Takahashi M, Jo S, Nishimura Y, Masunaga S, Fukuda M, Imura H. Radiofrequency thermotherapy for malignant liver tumors. Cancer. 1990;65(8):1730–1736. [PubMed]
23. Nau WH, Roselli RJ, Milam DF. Measurement of thermal effects on the optical properties of prostate tissue at wavelengths of 1064 and 633 nm. Lasers Surg Med. 1999;24(1):38–47. [PubMed]
24. Oden JT, Diller KR, Bajaj C, Browne JC, Hazle J, Babuška I, Bass J, Demkowicz L, Feng Y, Fuentes D, et al. Development of a computational paradigm for laser treatment of cancer. Lect Notes Comput Sci. 2006;3993:530. [PubMed]
25. Oden JT, Diller KR, Bajaj C, Browne JC, Hazle J, Babuška I, Bass J, Demkowicz L, Feng Y, Fuentes D, Prudhomme S, Rylander MN, Stafford RJ, Zhang Y. Dynamic data-driven finite element models for laser treatment of prostate cancer. Numer Meth PDE. 2007;23:904–922. [PMC free article] [PubMed]
26. Pennes HH. Analysis of tissue and arterial blood temperatures in the resting forearm. J Appl Physiol. 1948;1:93–122. [PubMed]
27. Robinson DS, Parel JM, Denham DB, Gonzalez-Cirre X, Manns F, Milne PJ, Schachner RD, Herron AJ, Comander J, Hauptmann G. Interstitial laser hyperthermia model development for minimally invasive therapy of breast carcinoma. J Am Coll Surg. 1998;186(3):284–292. [PubMed]
28. Roemer RB. Applications of bioheat transfer simulations in hyperthermia. Cancer Res. 1984;44(10):4788–4798. [PubMed]
29. Roemer RB, Dutton AW. A generic tissue convective energy balance equation: Part 1 – Theory and Derivation. J Biomech Engrg. 1998;120:395–404. [PubMed]
30. Rylander MN. PhD Thesis. The University of Texas; Austin: 2005. Design of Hyperthermia Protocols for Inducing Cardiac Protection and Tumor Destruction by Controlling Heat Shock Protein Expression.
31. Rylander MN, Feng Y, Bass J, Diller K. Hsp expression and damage optimization algorithm for prostate cancer therapy design. Lasers Surg Med. 2007;39:731–746. [PubMed]
32. Rylander MN, Feng Y, Bass J, Diller KR. Thermally induced injury and heat-shock protein expression in cells and tissues. Ann NY Acad Sci. 2005;1066:222–242. [PubMed]
33. Rylander MN, Feng Y, Zhang J, Bass J, Stafford RJ, Hazle J, Diller K. Optimizing hsp expression in prostate cancer laser therapy through predictive computational models. J Biomed Optics. 2006;11(4):04111:3-1–04111:3-16. [PubMed]
34. Showalter RE. Monotone Operators in Banach Space and Nonlinear Partial Differential Equations. American Mathematical Society; 1997.
35. Weinbaum S, Jiji LM, Lemons DE. Theory and experiment for the effect of vascular microstructure on surface tissue heat transfer. Part I. Anatomical foundation and model conceptualization. ASME J Biomech Engrg. 1984;106:321. [PubMed]
36. Weinbaum S, Xu LX, Zhu L, Ekpene A. A new fundamental bioheat equation for muscle tissue – Part I: Blood perfusion term. J Biomech Engrg. 1997;119:278. [PubMed]
37. Welch AJ, van Gemert MJC. Optical-Thermal Response of Laser-Irradiated Tissue. Plenum Press; New York: 1995.
38. White JA, Dutton AW, Schmidt JA, Roemer RB. An accurate, convective energy equation based automated meshing technique for analysis of blood vessels and tissues. Int J Hyperther. 2000;16(2):145–158. [PubMed]
39. Wissler EH. Pennes’ 1948 paper revisited. J Appl Physiol. 1998;85(1):35–41. [PubMed]
40. Wulff W. The energy conservation equation for living tissue. IEEE Trans Biomed Engrg. 1974;BME-21:494–495.
41. Wust P, Hildebrandt B, Sreenivasa G, Rau B, Gellermann J, Riess H, Felix R, Schlag PM. Hyperthermia in combined treatment of cancer. Lancet Oncol. 2002;3(8):487–497. [PubMed]
42. Yang R, Sanghvi NT, Rescorla FJ, Kopecky KK, Grosfeld JL. Liver cancer ablation with extracorporeal high-intensity focused ultrasound. Eur Urol. 1993;23(1):17–22. [PubMed]
43. Zhu L, Xu LX, He Q, Weinbaum S. A new fundamental bioheat equation for muscle tissue – Part II: Temperature of SAV vessels. J Biomech Engrg. 2002;124:121. [PubMed]
44. Zhu M, Weinbaum S, Jiji LM, Lemons DE. On the generalization of the Weinbaum–Jiji bioheat equation to microvessels of unequal size; the relation between the near field and local average tissue temperatures. J Biomech Engrg. 1988;110(1):74–81. [PubMed]