|Home | About | Journals | Submit | Contact Us | Français|
Laser surgery, or laser-induced thermal therapy, is a minimally invasive alternative or adjuvant to surgical resection in treating tumors embedded in vital organs with poorly defined boundaries. Its use, however, is limited due to the lack of precise control of heating and slow rate of thermal diffusion in the tissue. Nanoparticles, such as nanoshells, can act as intense heat absorbers when they are injected into tumors. These nanoshells can enhance thermal energy deposition into target regions to improve the ability for destroying larger cancerous tissue volumes with lower thermal doses. The goal of this paper is to present an integrated computer model using a so-called nested-block optimization algorithm to simulate laser surgery and provide transient temperature field predictions. In particular, this algorithm aims to capture changes in optical and thermal properties due to nanoshell inclusion and tissue property variation during laser surgery. Numerical results show that this model is able to characterize variation of tissue properties for laser surgical procedures and predict transient temperature fields comparable to those measured by in vivo magnetic resonance temperature imaging techniques. Note that the computational approach presented in the study is quite general and can be applied to other types of nanoparticle inclusions.
The latest statistics shows that cancer remains one of the leading causes of death in the United States . However, advances in nanotechnology and its applications in biomedical science and engineering over the past two decades have enabled numerous innovative and more effective cancer diagnosis and treatment modalities [3, 7, 10]. Treatment by traditional surgical resection procedures are a surgeon's standard approach for removal of well-defined primary tumors in non-vital regions. This technique is extremely invasive and usually associated with high morbidity. In contrast, thermal therapies employing a variety of heat sources including laser, focused ultrasound, and microwaves have benefits over conventional cancer treatment alternatives [12, 13, 19, 25, 26]. These treatment therapies are minimally invasive and can provide an alternative option to treat solid tumors embedded in vital regions. Technological advancements, such as actively cooled applicators and high power diode lasers, have made laser-induced thermal therapy more efficient, economical, and safer than other thermal therapeutic modalities. Some advantages include that laser-induced thermal therapy can be used to treat tumors more rapidly than other modalities and have more control over perfusion effects. Additionally, lasers do not require a complicated setup that involves grounding pads and can be incorporated safely into any imaging environment, including MRI, with minimal induced image artifacts. The interaction between multiple probes for treating larger tumors can be synergistic and fully compatible with MRI for online monitoring.
On the other hand, nanoshell-mediated laser surgery is able to direct thermal energy into target regions delivered by optical fibers to provide a lethal dose of heat while minimizing damage to surrounding tissue . In particular, laser surgery promises effective treatment of small, poorly-defined metastases or other tumors embedded within vital regions. In this study, we consider a special class of nanoparticles known as nanoshells, which can act as intense infrared absorbers which increase the thermal deposition of laser energy into tumors. In particular, nanoshells provide a potential means to (a) enhance the delivery of laser-induced thermal energy via distributing the heat source from the fiber to the surrounding vasculature and/or, (b) 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. Typically nanoshells consist of a concentric spherical dielectric (silica) core and a thin metal coating (Au) shell. The diameters of nanoshells are usually in the 110 to 120 nm range and have been shown to be effective mediated agents to control the temperature field. Nanoshells possess a highly tunable plasmon resonance which determines the particle's scattering and absorbing properties. The plasmon resonance, one of the nanoshell's optical properties, can be tuned across a broad range of the light spectra from the ultra-violet to the infrared by controlling the ratio between the radius of the core and the thickness of the shell layer [1, 11]. When nanoshells are injected to the target region, laser-induced thermal energy can be delivered to specified locations and greatly enhance heat absorption in tumor regions due to the change of optical properties in the tissue .
To design optimal nanoshell-mediated laser surgical protocols, it is crucial to accurately characterize the optical, thermal, and biological response of tissues to therapies [20, 21]. The major challenge is that these properties can be difficult to measure and vary over time during the treatment due to biological alteration in tissues. The goal of this paper is to present a novel nested-block optimization algorithm for nonlinear transient bioheat transfer model to simulate laser surgery in the presence of nanoshells and with consideration of dynamic changes in optical and thermal properties of the tissue due to biological alteration. In this study, a three-dimensional finite element nonlinear transient bioheat transfer model with input of laser-tissue interaction calculation from Monte Carlo fluence model is constructed. Numerical results show that this model can reliably characterize changes in tissue properties and accurately predict temperature fields comparable to those measured by in vivo magnetic resonance temperature imaging (MRTI) technique. Although the validation experiments are conducted for treating prostate tumors inoculated on SCID (severely compromised immuno-deficient) mice, the computational approach presented in the study is quite general and can be applied to other types of cancer treatment. Similar optimization strategies are also used in a dynamic data-driven framework for real-time surgical control .
Human prostate cancer cells (PC-3) were cultured with HAM's F12 medium with 10% FBS and 5% penicillin-streptomycin. Mice were anesthetized with isoflurane. PC-3 cells were then inoculated in the backs of 4–6 week old SCID mice with an injection volume of 0.2 ml of 5 × 106 PC-3 cells in Matrigel for each mouse. Prostate tumors were grown to a tumor burden of approximately 1.0 cc (equivalently a spherical tumor with radius of 6.2 mm). Figure 1a depicts a typical PC-3 tumor on a SCID mouse prepared for laser irradiation.
Similar to previous work , nanoshells were employed in conjunction with extracorporeal laser irradiation to enhance thermal deposition. Nanoshells were injected into the mouse tail vein 24 h prior to laser irradiation to enable adequate accumulation in the tumor volume. Nanoshells (furnished by Nanospectra Biosciences Inc. Houston, TX) are composed of a silica core (with a diameter of 110 nm) and an outer gold shell (thickness of 15 nm). The shell geometry was optimized to enable maximum absorption at wavelength 808 nm to enhance thermal deposition. Figure 1b illustrates the distribution and geometry of the nanoshells employed in laser therapy experiments.
By tuning the ratio of the core diameter and shell thickness, nanoshells can be made to absorb or scatter light at a desired wavelength across visible and near-infrared (NIR) wavelengths. This optical tunability permits optimal design of laser surgical protocols with a peak optical absorption in the NIR region. For instance, laser surgery for deeper tissue requires light in the NIR region where tissue has the highest transmissivity. The ability of gold-nanoshells to convert strongly absorbed light into localized heat can be exploited for the targeted laser therapy of cancer. Thus, effective targeting of nanoshell bioconjugates specifically to cancer cells combined with the high absorption cross-section of the nanoshells in the laser excitation region, generates increased temperatures sufficient to produce irreversible cell and tissue damage to subcutaneous tumors while keeping laser energy at a lower level so that cells outside of the target region are minimally damaged.
During the laser surgical experiments, the temperature distribution was measured by in vivo MRTI with the proton-resonance frequency-shift method [8, 17]. All experiments were performed at the University of Texas M.D. Anderson Cancer Center in Houston, Texas, on a 1.5-T MR scanner (Excite HD, GEHT, Waukesha, WI) equipped with high-performance gradients (23 mT m–1 maximum amplitude and 120 T m–1 s–1 maximum slew rate) and fast receiver hardware (bandwidth, ±500 MHz). Mice were imaged with a three-inch receive-only surface coil specially designed for small animal imaging. T1- and T2-weighted images were used to plan and localize the treatment by verifying the position of the laser fiber relative to the imaged region prior to irradiation. The skin over the tumor was swabbed with PEG diacrylate (Mx 600, Sartomer, West Chester, PA) as an index-matching agent. MRTI was performed by using a complex phase-difference technique with a fast, two-dimensional RF-spoiled gradient-recalled echo sequence (TR/TE = 49.5/20 ms, flip angle = 30°, bandwidth = 9.62 kHz). To achieve a five-second per image scan rate, a rectangular field of view (6 × 4 cm2) with 256 × 86 encoding matrix was used. Also the reduced bandwidth was employed to minimize gradient heating limitations at this small field and enhance the signal-to-noise ratio. The acquired voxel size was 1.07 × 0.82 × 3 mm3. The change in temperature from baseline after a number of images was extrapolated from the complex-valued MRTI data by using the temperature dependence of the proton resonance frequency shift and a temperature sensitivity with coefficient of –0.0097 ppm/°C . The temperature resolution for MRTI measured temperature difference is less than 1°C. Figure 2 shows typical MRI and MRTI images along with a quantified image that combines the two.
Simulating laser surgery and making reliable predictions of the temperature field requires two major modeling components: a bioheat transfer model for the tissue and a laser source term that characterizes thermal energy deposited into the tissue. Since the optical properties of nanoshells can be designed by adjusting the ratio between the diameter of the silica core and the thickness of the gold shell, tissue properties such as the absorption and scattering coefficients can be controlled then with inclusion of nanoshells. In this section, we discuss the mathematical and computational models for bioheat transfer and laser-tissue interaction.
The mathematical representation of the temperature distribution in the tissue incorporates both the Pennes bioheat equation for the thermal effects of local blood perfusion and an expression for laser energy as a thermal source. We consider the case of external heating provided by a laser source on the surface of the tumor.
The thermal conductivity, k [J s–1 m–1 K–1], and blood perfusivity, ω [kg s–1 m–3], are assumed to be bounded functions of the temperature field, T = T(x, t), where
and k0[J s–1m–1 K–1], k1 [J s–1m–1K–1], k2[K–1], k3[K], ω0[kg s–1 m–3], ω1[kg s–1 m–3], ω2[K–1], and ω3[K] are parameters of the diffusivity and perfusivity coefficient functions defined above. The specific heat of the tissue and blood are given by cp and cb respectively, Ta is the arterial temperature, and ρ is the density of the tissue. On the Cauchy boundary,
The coefficient of cooling is denoted h. T∞ denotes the ambient temperature. is the prescribed heat flux on the Neumann boundary.
The temperature field is propagated forward in time from a given initial condition, T0.
Laser-tissue interaction is a complex phenomenon and usually divided into optical and thermal responses [15, 24]. Typically, laser light is absorbed by tissue and converted into heat. A commonly used model that characterizes the absorbed heat can be represented by the rate of heat generation Q(x, t) that is defined as
and Φ(x, t) is the fluence that defines the amount of energy, in the form of photons, passing through a unit area at a point in space per unit time. P(t) is the laser power at time t. The parameters μa, μs are absorption and scattering coefficients that relate to laser wavelength. They can be considered as the probability of absorption and scattering, respectively, of photons in tissues. The constant g is the anisotropic factor, and x0 is the position of the laser source.
Derived from the light diffusion approximation for monoenergetic neutral particles by solving the transport equation , (2) is valid when the scattering coefficient is much larger than the absorption coefficient (), which is the case in most of the soft tissues. The basic assumption in our case is that the light is emitted from a single point in a diffuse and isotropic manner. In addition, this model only accounts for scattered light, i.e., the primary light is ignored. This assumption is reasonable if the region of interest is far from the source. For external heating, these assumptions may be violated. In our case, the light from the laser beam was collimated as it enters the tissue and is incident against the skin with a flat beam spot between 0.5 mm and 1.0 cm in diameter.
In this study, we use both the classical isotropic approach, (2), and the Monte Carlo method to solve for the fluence (see e.g. ). Generally speaking, the Monte Carlo method is considered to be more accurate but computationally more expensive than the classical approach. We use the Monte Carlo method here to verify the accuracy of the classical isotropic approach. The two major assumptions that are used in the Monte Carlo model are: (a) the photons will be distributed in a cylindrically symmetric manner with the initial direction of the beam as the axial direction, and (b) the probability distributions being used for the scattering angles and mean free path length are known. To compute the fluence, we consider the following three quantities of interest as usual.
All three quantities depend on light wavelength, temperature, and space, due to heterogeneity of the tissue. However, the space dependency is usually ignored. Thus, μa, μs, and g are typically given as functions of light wavelength and temperature only. It has been found that values of these three quantities remain fairly constant for temperatures below 70°C .
The Monte Carlo simulation for laser irradiation starts with following the random paths of many individual photons by sampling probability distributions. The main idea is described as follows. The algorithm starts by initiating a photon with a weight of W = 1 and a position and direction in space. A random number between 0 and 1 is then generated and used as a value of a mean free path length by correlating it to the associated probability distribution. The photon is then moved with this length in its current direction. Then, the algorithm checks if the photon is still in the tissue. If it is not, the internal reflection for this particular photon is done. Otherwise, a portion of the photon, W · μa/μt, is assumed to be absorbed at that location, the weight is updated and then two more random numbers are generated to obtain the new azimuthal angle and the deflection angle by again correlating the random numbers to the respective probability distribution. This process is repeated until the photon's weight has reached some cutoff weight. To conserve energy, a roulette scheme is employed here to decide if the photon should be discarded or left in with an update of its weight.
The output of Monte Carlo algorithm is the number of photons absorbed in each grid cell descretized in an axi-symmetric r–z plane. Then, the fluence can be obtained by
where grid[z, r] is the number of photons absorbed in the grid cell [z, r], M is the total number of photons put into the simulation, V[z, r] is the volume represented by the grid cell [z, r], and μa is the absorption coefficient. The grid size can be refined according to desired tolerance. Having computed Φ[z, r] at each discrete grid cell, it will be projected to the three-dimensional finite element mesh and used to calculate the heat generation term Q(x, t) = μaΦ(x, t) in (1).
In this section, we present a nested-optimization algorithm that applies to the Pennes bioheat transfer model . The goal is to capture dynamic changes in tissue properties due to biological alteration as well as nanoshell inclusion. In order to achieve this goal, the objective function, (5), for the calibration problem is defined as the difference in space-time norm between the computed temperature field, T(x, t), and the temperature field obtained from the MRTI experiments, Texp(x, t), integrated over the entire biological domain Ω and time duration of interest [0, τ]. The computed temperature field T(x, t) is an implicit function of the bioheat transfer model parameters β, a vector that consists of the location of the laser probe, absorption and scattering coefficients, and parameters in conductivity and perfusion functions.
The main question to be addressed in this section is how to approximate the full optimization problem efficiently so that the prediction can be made prior to the next data arrival. In other words, the speed of computation needs to be faster than the rate of data acquisition, at least faster than the time duration of the operation. The optimization problem can be formally stated as:
Find β* such that Q(T(β), β) is minimized, where
In order to speed up solution time for the approximation, (3) is replaced by a sequence of smaller problems with less historical data in the time dimension.
Find β*i such that Qi(T(β),β) is minimized, where
using a smaller time window , i = 1,2,3,.... Figure 3 illustrates the nested-block optimization process: Calibration → Prediction → Validation/Calibration → Prediction, as in vivo MRTI measurement data continue to be imported into the model.
Based on the experiment, the temperature field is calibrated with respect to MRTI thermal images of heating of a mouse tumor treated with nanoshells. A nonlinear form of Pennes bioheat transfer model, (1), is calibrated to predict the heat transfer in tissues. There is a trade-off between accuracy and efficiency. An efficient optimization process will allow faster solution time for the calibration so that it can be used in real-time surgical control. However, time constraints of real-time computation do not permit the computation of the Hessian of the objective function, which would allow a more accurate measure of rate change with respect to each model parameter as a solution is approached. Thus, we compute the gradient of the objective function using a limited-memory variable metric by a quasi-Newton optimization method . As shown in , the gradient vector of the quantity of interest can be written as
where  denotes  dx dt and p(x, t) is the solution to the adjoint problem
with the boundary conditions
and the terminal condition
The major computational modules and data flow for the integrated nonlinear transient bioheat transfer computational model with laser heating are shown in Fig. 4.
The nonlinear transient bioheat transfer model is implemented using a finite element discretization in space and finite difference in time using the Crank-Nicolson scheme.
The laser source model was implemented using both the isotropic analytical and the Monte Carlo models as described previously. For the Monte Carlo calculation, results were obtained by averaging over 10 runs with 1,000,000 photons in each run. Parameters used were μa = 2.15 cm–1 and μs = 14.2 cm–1. The grid size was chosen Δr = 0.015 cm–1 and Δz = 0.005 cm–1. The cumulative distribution function for the mean free path length is defined as
where μt = μa + μs. We assume that the probability distribution of the azimuthal angle is uniform. The probability distribution of the cosine of the deflection angle is defined as
which is known as the Heyney-Greenstein scattering function, where the anisotropic factor g is set to 0.9 (a typical value for soft tissues).
Both bioheat transfer and optimization modules are implemented in parallel using an h–p finite element method built upon an in-house code and can be executed efficiently on a multi-processor machine. The framework of this adaptive finite element code has been developed and thoroughly tested over last two decades . Computations were carried out using Linux clusters hosted at the Texas Advanced Computing Center at Austin. Each computing node has two Xeon Intel Duo-core 64-bit/2.66GHz processors. Using up to 50 processors, we are able to make real-time predictions two minutes ahead, which only takes 10 s computing time.
Data acquired for calibration are shown in Fig. 2. Figure 2a is a 49 × 56 pixel MRI image of a mouse tumor. The field of view is 4 × 6 cm2 and the thickness associated with the image is 3 mm. An external heat source was applied to a mouse tumor. Sixty thermal images were acquired at every 5 s. A single time instance of the thermal images is shown in Fig. 2c. An overlay of the thermal image onto the geometry is shown in Fig. 2b.
In order to compare computational results with the experimental data, the in vivo temperature measurement (MRTI) data are interpolated onto a finite element mesh of the biological domain generated from the MRI data. In Fig. 5, the results at time t = 300 s are shown. The initial mesh consists of 6,076 linear elements with a total of 8,000 degrees of freedom. The MRTI data was filtered and then projected onto the finite element mesh, which is shown in Fig. 5a. Figure 5b compares a cut-line through the mid-plane of the unfiltered thermal image to the filtered thermal image. The Crank-Nicolson scheme was used to propagate the bioheat transfer equations forward in time. Figure 5c and e show the thermal contour of the FEM prediction corresponding to the isotropic laser source term, (2), and the Monte Carlo input respectively. Parameters used in the simulation are given in Tables 1 and and2.2. The initial values at the beginning of the calibration process are taken from .
Figure 5d and f illustrate a cut-line comparison of the MRTI data to the FEM prediction using the isotropic laser source term and the Monte Carlo source term. The FEM model with uncalibrated parameters seems to produce inaccurate results compared with MRTI data, which strongly indicates that the tissue properties are changed due to the introduction of nanoshells.
An early study  shows that the model of bioheat transfer with a laser heating source term is anticipated to be relatively less sensitive to the attenuation coefficients of the laser source term as compared to the thermal conductivity and blood perfusivity coefficients. Thus, more changes in the blood perfusivity and tissue conductivity are expected than in absorption and scattering coefficients in the calibration process.
Using nested-block optimization algorithm, calibrations were performed using different sets of thermal images in the objective function with number of images N = 10, 20, 30, 40, 50, 60. The variation in the calibrated coefficients is shown in Fig. 6a–d. It is observed that the perfusion coefficients and thermal conductivity parameters are the most sensitive to the calibration. Figure 6e and f shows the sensitivity of the objective function and its gradient to the number of thermal images used in the optimizations. It was observed that the objective function is rather insensitive to the number of thermal images used in the optimizations of N > 30. For time critical applications, this implies that the calibration process may be just as accurate using only half the full data set. For comparison at a different time instance, thermal contour and a cut-line of the FEM predictions using N = 20 and 60 is given in Fig. 7 at time t = 235 s.
The extent of thermal damage associated with traditional laser therapies is limited by thermal diffusion in the tissue. Most conventional laser therapies, other than photodynamic therapy, do not permit selective destruction of tumor tissue while preserving the surrounding healthy tissue. In order to circumvent the limitations of traditional laser therapies employed in prostate cancer treatment, nanoshells were employed to enhance the thermal deposition and selectivity of thermal damage.
In this study, a new computational methodology for solving the Pennes bioheat transfer equation which incorporates a nested-block optimization algorithm is proposed. This methodology is used to simulate the transient temperature field during laser surgery on a prostate tumor inoculated on the back of SCID mice. The numerical simulations for this test case are in good agreement with the MRTI data from the surgical procedure but required several iterations through the calibration process. For the single test case presented herein the calibration process stabilized after 30 time steps (150 seconds into the surgery).
A flexible Monte Carlo based formulation for the fluence in the laser source term is incorporated in the Pennes bioheat equation, which is solved by a finite element method. This formulation is considerably more general than the classical approach allowing for non-homogeneous and anisotropic tissues, and extends the range of validity to include more complex tumor geometry. Preliminary results using this more general laser model show only modest variation (lest than 5%) in the temperature field when compared to results produced using a classical formulation for the fluence. We anticipate the new formulation will be a more significant factor in the analysis as we investigate more complex tumor geometry from the perspective of laser source modeling. Additional validation tests involving more specimens with different laser treatments protocols are underway to further refine this methodology.
The computer model developed in this study for laser-induced thermal therapy is capable of predicting the temperature due to laser heating in a prostate tumor with nanoshell inclusion by incorporating changes in the absorption and scattering effects. Using this model, optimal hyperthermia protocols can be developed. With maturity of nanotechnology such as nanoshell injection and a reliable computational model, the outcome of laser surgery for cancer treatment will become more controllable and predictable. MRTI provides an in vivo measurement tool for mapping the spatiotemporal temperature elevation in the tumor during and following laser therapy. This technique for temperature measurement has resolution of less than 1°C. With the nested-block optimization algorithm presented in this study, tissue properties can be accurately characterized during the course of a surgical operation and resulting in more reliable prediction of temperature field and treatment outcomes.
The support of this work by the National Science Foundation under grant CNS-0540033 and by National Institutes of Health under Small Animal Imaging Facility core grant CA 16672 are gratefully acknowledged.
Yusheng Feng, Computational Bioengineering and Nanotechnology Lab, Department of Mechanical Engineering, The University of Texas at San Antonio, San Antonio, TX 78249, USA ; Email: email@example.com.
David Fuentes, Institute Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA.
Andrea Hawkins, Institute Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA.
Jon Bass, Institute Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA.
Marissa Nichole Rylander, Department of Mechanical Engineering and School of Biomedical Engineering and Sciences, Virginia Tech, Blacksburg, VA 24061, USA.
Andrew Elliott, Department of Imaging Physics, The University of Texas M.D. Anderson Cancer Center, Houston, TX 77030, USA.
Anil Shetty, Department of Imaging Physics, The University of Texas M.D. Anderson Cancer Center, Houston, TX 77030, USA.
R. Jason Stafford, Department of Imaging Physics, The University of Texas M.D. Anderson Cancer Center, Houston, TX 77030, USA.
J. Tinsley Oden, Institute Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA.