Nanoshell-mediated laser surgery simulation for prostate cancer treatment

2.1 Tumor preparation and inoculation

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 × 10⁶ 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.

2.2 External heating and nanoshell inclusion

Similar to previous work [9], 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.

2.3 Temperature measurement

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. T₁- and T₂-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 (M_x 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 cm²) 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 mm³. 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 [8]. 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.

3.1 Bioheat transfer model

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.

Let Ω be a bounded domain (tumor region) in ${\mathcal{R}}^3$ with $\Gamma =\partial {\Omega }_{\mathrm{C}}\bigcup \partial {\Omega }_{\mathrm{N}}$ denoting a Lipschitz continuous boundary. Equation (1) is based on Pennes bioheat equation [18] with temperature dependent coefficients:

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 k₀[J s^–1m^–1 K^–1], k₁ [J s^–1m^–1K^–1], k₂[K^–1], k₃[K], ω₀[kg s^–1 m^–3], ω₁[kg s^–1 m^–3], ω₂[K^–1], and ω₃[K] are parameters of the diffusivity and perfusivity coefficient functions defined above. The specific heat of the tissue and blood are given by c_p and c_b respectively, T_a 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. $\mathcal{G}$ is the prescribed heat flux on the Neumann boundary.

The temperature field is propagated forward in time from a given initial condition, T₀.

3.2 Laser-heating model

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 x₀ is the position of the laser source.

Derived from the light diffusion approximation for monoenergetic neutral particles by solving the transport equation [23], (2) is valid when the scattering coefficient is much larger than the absorption coefficient (${\mu }_s\gg {\mu }_a$), 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. [23]). 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.

1. the absorption coefficient, μ_a, the average number of photons absorbed per unit length,

2. the scattering coefficient, μ_s, the average number of photons scattered per unit length, and

3. the anisotropic factor, g, the expected value of the cosine of the deflection angle.

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

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

3.3 Nested-block optimization algorithm

In this section, we present a nested-optimization algorithm that applies to the Pennes bioheat transfer model [18]. 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, T_exp(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 Q_i(T(β),β) is minimized, where

using a smaller time window $[0,{\tau }_i]\subset [0,{\tau }_{i+1}]\subset [0,\tau ]$, 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 [2]. As shown in [16], the gradient vector of the quantity of interest can be written as

where [] denotes ${\int }_0^{\tau }{\int }_{\Omega }$ [] dx dt and p(x, t) is the solution to the adjoint problem

with the boundary conditions

and the terminal condition

3.4 Computational implementation

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