|Home | About | Journals | Submit | Contact Us | Français|
The treatment times of laser induced thermal therapies (LITT) guided by computational prediction are determined by the convergence behavior of PDE constrained optimization problems. In this work, we investigate the convergence behavior of a bioheat transfer constrained calibration problem to assess the feasibility of applying to real-time patient specific data. The calibration techniques utilize multi-planar thermal images obtained from the non-destructive in vivo heating of canine prostate. The calibration techniques attempt to adaptively recover the bio-thermal heterogeneities within the tissue on a patient specific level and results in a formidable PDE constrained optimization problem to be solved in real time. A comprehensive calibration study is performed with both homogeneous and spatially heterogeneous bio-thermal model parameters with and without constitutive nonlinearities. Initial results presented here indicate that the calibration problems involving the inverse solution of thousands of model parameters can converge to a solution within three minutes and decrease the norm of the difference between computational prediction and the measured temperature values to a patient specific regime.
The feasibility of developing an adaptive feedback control system which uses dynamic multi-planar MRI temperature measurements to adaptively conform delivery of therapy to the prescribed treatment plan in real-time during MR-guided laser induced thermal therapy has been demonstrated . The cyberinfrastructure ,  inherent to the control system relies critically on the precise real-time orchestration of large-scale parallel computing, high-speed data transfer, a diode laser, dynamic imaging, visualizations, inverse-analysis algorithms, registration, and mesh generation. We demonstrated that this integrated technology has significant potential to facilitate a reliable minimally invasive treatment modality that delivers a precise, predictable and controllable thermal dose prescribed by oncologists and surgeons. However, MR guided LITT (MRgLITT) has just recently entered into patient use  and substantial translational research and validation is needed to fully realize the potential of this technology ,  within a clinical setting. The natural progression of the computer driven MRgLITT technology will begin with prospective pre-treatment planning. Future innovations on the delivery side will likely involve combining robotic manipulation of fiber location within the applicator as well as multiple treatment applicators firing simultaneously.
Currently, the workflow of our treatment protocols are divided into data acquisition, computational, and delivery phases. The present work focuses on the computational aspects of model calibration. A schematic of the procedure setup and workflow is provided in Figure 1 and Figure 2, respectively. Heterogeneous model calibration involving thousands of model parameters have been shown to deliver model predictions of unprecedented accuracy . The goal of this study is to evaluate the feasibility of real-time patient specific heterogeneous model calibration for both linear and nonlinear constitutive models. Such calibrations are critical for maintaining the predictive power of the simulation during therapy and therefore in maximizing the efficiency of the therapy control loop. We briefly review the governing PDE constrained optimization equations, present our results, and discuss the higher level implications for real-time control of MR-guided LITT as well as milestones within the translational research.
The problem of bioheat transfer model calibration is to determine the set of thermal parameters that minimize the L2(0, T; L2(Ω)) norm of the difference between the predicted temperature field, u(x, t) and the temperature field observed in vivo thermal images of the experiment, uMRTI(x, t). In our experiments, the MR temperature imaging (MRTI) data is acquired using a 2D multi-slice temperature sensitive echo planar imaging sequence collecting five planes of temperature sensitive images every five seconds  on a 1.5 T MRI scanner. The cost function is then
dx = dx1dx2dx3 being a volume element and the time interval of interest is denote ΔT. The simulation of the time evolution of temperature field within the biological domain is constrained by the classical Pennes model  of bioheat transfer with a isotropic laser heat source.
Pennes model has been shown to provide very accurate prediction of bioheat transfer , ,  and is used as the basis of the finite element prediction. The full initial boundary value model is defined by the following system:
The initial temperature field, u(x, 0) = u0, is taken as the measured baseline body temperature. The density of the continuum is denoted ρ and the specific heat of blood is denoted . On the Cauchy boundary, ΩC, h is the coefficient of cooling and u∞ is the ambient temperature. The prescribed heat flux on the Neumann boundary, ΩN, is denoted . The optical-thermal response to the laser source, Qlaser(x, t), is modeled as the classical spherically symmetric isotropic solution to the transport equation of light within a laser-irradiated tissue . P(t) is the laser power as a function of time, μa and μs are laser coefficients related to laser wavelength and give probability of absorption and scattering of photons, respectively. The anisotropic factor is denoted g and x0 denotes the position of laser photon source. The scalar-valued coefficient of thermal conductivity is modeled with a nonlinear temperature relation
where , k3 [K] . Perfusion is modeled with a nonlinear dependence on temperature, Figure 3:
where , ωNI [K], ωID [K] . The assumed perfusion model is an extension of the constitutive model used in  and attempts to recover the observed physiological effects of perfused tissue. During a thermal therapy, the applied heat begins to dilate the vasculature at a temperature of ωNI and the normal value of the perfusion, ωN, increases to a state of hyper perfusion, ωI. Beyond a critical threshold temperature, ωID, the vasculature is damaged and a reduced amount of perfusion is seen, denoted ωD. In addition to modeling the vasculature breakdown above a temperature threshold, the perfusion model attempts to capture the expected hyper-perfusion under hyperthermia conditions. The linear components of the perfusion, ω0(x), and thermal conductivity, k0(x), are allowed to vary spatially within a local region of interest, r = 1cm, around the laser source.
Constitutive model data used is summarized in Table I. For the data in this study, the heating is localized to a ≈1cm diameter region about the diffusing interstitial fiber tip and is not significantly influenced by the boundary during delivery. Hence, zero flux boundary conditions are used, = 0 and ΩC = . The perfusion and thermal conductivities, highlighted in (2), are recovered from the calibration computation. The linear heterogeneous model and the nonlinear constitutive models are included in the parameter space explored by the optimization scheme.
We employ a limited-memory variable metric that arises in a quasi-Newton optimization method  to drive the calibration problem. Using an adjoint method, the gradient of the objective function (1) with u the solution of (2) can be written as 
Here p is the solution to the adjoint problem and 0(x), 1, 2, 3, 0(x), N, I, D, 2, NI, ID are model parameter test functions.
In vivo MR-guided LITT experiments were performed at The University of Texas M.D. Anderson Cancer Center in Houston, Texas. Handling of the canine was in accordance with an Institutional Animal Care and Use Committee approved protocol. General anesthesia was induced utilizing meditomidine (0.5 mg/kg, intramuscular) and 2% isoflurane was used to maintain general anesthesia throughout the duration of the experiment. The experimental configuration for the calibration study is illustrated conceptually in Figure 1. A FEM mesh of the full prostate anatomy for treatment planning, Figure 1, was created using pre-operative axial imaging data of the canine prostate and the neighboring anatomy. Two first order hexahedral meshes were used in the study. The lower resolution mesh consisted of 8817 elements and 9344 nodes and the higher resolution mesh consisted of 22376 elements and 24539 nodes. The meshes consisted of 8 and 12 elements across the diameter of the heating region, respectively. An example of the lower mesh resolution is provided in the cutplanes shown in Figure 4(a).
Axial and coronal planning images were acquired and used in conjunction with fiducials on a planning template to guide the position of the laser fiber before applying the power. All images were acquired on a clinical 1.5-T MR scanner (Excite HD, GEHT, Waukesha, WI) equipped with high-performance gradients (40 mT/m maximum amplitude and 150 T m−1 sec−1 maximum slew rate) and fast receiver hardware (bandwidth, ±500 MHz). A stainless steel stylet was used for inserting the laser catheter, 400 micron core diameter silica fiber in a water-cooled diffused tip catheter (Visualase Inc. Houston, TX). The location of the laser, in DICOM coordinates, was established from the intra-operative images and used in the calibration simulations. A region of the prostate of an anesthetized dog was heated with a nondestructive test pulse from an interstitial laser fiber (980nm; 5 Watts for 60 seconds) housed in an actively cooled applicator Δt. The diffusing interstitial laser source is expected to emit a uniform and homogeneous fluence.
Real-time multi-planar MRTI based on the temperature dependence of the proton resonance frequency (PRF)  monitored the heating and cooling phases of the laser test pulse. The on/off times of the laser were known and controlled by the software so that a priori knowledge of beam on and off time was known. The PRF technique uses a complex phase subtraction technique; each consecutive phase image is subtracted from a baseline image that was acquired prior to heating . These phase difference images are proportional to temperature difference, Δu. To convert phase to temperature, other factors including echo time (15ms), water proton resonance frequency (63.87MHz), and the temperature sensitivity of the water proton (−0.0097 ppm/°C ) are required. The image acquisition time of the thermal imaging data in this study is 5 planes of complex image data every 5 seconds. In addition to a spatial median-deriche filtering pipeline, a space-time filter was also applied to the thermal imaging data; if the thermal data at a pixel changes by more than 11°C it is considered noise and filtered. The 11°C threshold is based on rates of heating observed in previous in vivo experiments in conjunction with the a priori knowledge that we are acquiring via a low power test pulse.
A calibration study to refine simulation parameters was performed varying the amount of thermal imaging data used in the objective function (1) as a function of the heating pulse time Δt 60 seconds. The performance of the metric chosen for the optimization process was studied in terms of the time window of data used and the effect on the final value of the space time norm between the data and the model prediction over the time window ΔT = (0, 2.0Δt). The metrics considered were the space time norm between the data and the model over 5 time intervals of interest: ΔT = (0, 0.5Δt), (0, 1.0Δt), (0, 1.5Δt), (0, 2.0Δt), and (1.0Δt, 2.0Δt). The time intervals ΔT = (0, 0.5Δt) and (0, 1.0Δt) were chosen to study the effect of using only the heating phase for the calibration calculations. The time intervals ΔT = (0, 1.5Δt) and (0, 2.0Δt) captures the heating and cooling of the in vivo tissue. The time interval ΔT = (1.0Δt, 0, 2.0Δt) attempts to eliminate the laser source and isolate the tissue specific properties from the calibration problem and account only for the cooling of the tissue. Moreover, the behavior of the calibration problem was studied using homogeneous model coefficients versus heterogeneous model coefficients with and without the nonlinearities of the constitutive relationships for the perfusion and thermal conductivity.
A post nonlinear, heterogeneous calibration comparison between the Pennes model and thermal imaging data projected onto the finite element mesh is given in Figures 4 and and5.5. Figure 4 compares a cutplane that intersects the expected plane of highest heating pre- and post-calibration. The heterogeneous tissue property recovery allows for the non-isotropic heating contours observed in the Pennes model prediction, Figure 4(a). The relative position of the temperature profiles and temperature history, Figure 5, are shown in Figure 4(a). The temperature-time history during the cooling agrees with the exponential decay of the temperature expected from a classical separation of variables solution to the linear heat equation. The calibrated model prediction shows good agreement with the thermal images. The main source of disagreement in the pointwise error, Figure 4(a) and Figure 5(a), can be attributed to noise in the thermal imaging data away from the laser probe. However, the average standard deviation of pre-heating thermal images measured in a 5 x 5 x 5 pixel ROI within the contra-lateral lobe of the prostate was 1.53°C
The benefit of utilizing spatially heterogeneous techniques for model calibration is shown in Figures 6 and and7.7. Corresponding cutplanes through the uncalibrated model prediction with homogeneous coefficients is compared to the calibrated model prediction with heterogeneous coefficients, Figure 6(b). The corresponding thermal imaging data is provided for a reference, Figure 6(a). The model prediction with homogeneous coefficient is seen to create spherical isotherms; this is expected from an isotopic source. Model predictions with heterogeneous coefficients are seen able to recapitulate the structure of the isotherms in the thermal imaging data. This phenomenon is further illustrated in the temperature profiles provided in Figure 7. Profiles through the homogeneous case are symmetric about the axis of heating as expected, but profiles through the heterogeneous case are asymmetric and better agree with the thermal imaging data.
A comprehensive summary of the calibration study is provided in Figure 8. The study was done at two discretizations of the geometry, 9344 dof and 24539 dof using linear homogeneous model coefficients, nonlinear homogeneous coefficients, heterogeneous linear coefficients, and heterogeneous nonlinear coefficients. The linear homogeneous, nonlinear homogeneous, heterogeneous linear, and heterogeneous nonlinear cases resulted in, 2/2 dof, 11/11 dof, 937/3723 dof, 946/3731 dof for the optimization problem; respectively, at the two model discretizations studied. The space-time norm over the entirety of the data was used as the basis for comparison, i.e. for ΔT = (0, 0.5Δt) the model was calibrated using the objective function and then propagated overp the entire time interval ΔT = (0, 2.0Δt) as a metric for comparison against other data time intervals. The initial pre-calibration value of the objective function is provided as a reference. Results indicate that using a data interval ΔT = (0, 1.5Δt) for the calibration problem provides the greatest objective function decrease for the least amount of computational work and is thus optimal in terms of computation efficiently. The benefit for using the heterogeneous calibration over the homogeneous calibration is clearly in the difference of the magnitude of the objective function (1), however, the benefit of the constitutive nonlinearities is not as dramatic.
The observed objective function convergence history for selected calibration problems is presented in Figure 9. The graph is intended to convey the general convergence behavior of the breadth of the calibration problems studied. The number function evaluations required for convergence of the calibration problem for homogeneous and heterogeneous models coefficients with and without constitutive nonlinearities is shown for the (0, 0.5Δt), (0, 1.0Δt), (0, 1.5Δt), (0, 2.0Δt), and (1.0Δt, 2.0Δt) time windows. The insert of Figure 9 shows the- convergence trend of the optimizer for the linear homogeneous case with heterogeneous model coefficients at the highest resolution considered in this study. The two significant decreases in the objective function seen in Figure 9 were typically associated first with a large change in the thermal conductivity followed by a change in the perfusion. On average, the optimizer is seen to have converged by twenty function-gradient computations. As observed in the insert of Figure 9, the calibration problem appears to converge well before the typical relative and absolute value optimization convergence metrics are reached. The metric used to report the convergence in Figure 9 was weakened to a form (3) that retrospectively looks at the history and uses the relative difference between the absolute minimum-over converged function value obtained, Qminimum, and the starting value, Qinitial.
For a particular problem, when Qconverged satisfies (3), convergence is reached. Results provide an estimate that the lower bound wall clock time needed for a real-time patient-specific calibration computation is 3 minutes (for this particular study); 90 seconds worth of simulation for a function-gradient computation takes 9 seconds , multiplied by twenty objective function-gradient computations needed for convergence.
Calibration of the model based on accounting for spatial heterogeneities in the tissue produce much greater improvements in predicted temperature than upgrading terms in nonlinear constitutive equations assuming homogeneous tissue. This suggests an effective constitutive alternative to modeling the nonlinear bioheat transfer observed in soft tissue. In light of the fact that this work was limited to homogeneous constitutive nonlinearities, the benefit of allowing for a spatially varying constitutive nonlinearity in the calibration is not clear. Linear model heterogeneities provide a dramatic increase in the model prediction. The possible relative increase in model accuracy obtained by adding spatially varying constitutive nonlinearities does not seem to be worth the cost of implementation. The tissue behaves as a heterogeneous but linearly conductive media.
The thermal source in this study was modeled as isotropic. As illustrated by the agreement in the model prediction, Figure 6, the modeling error in the laser source appears to have been compensated by the spatially varying thermal parameter field. Further work is needed to differentiate the extent of effects of the actual tissue heterogeneity, the photon distribution emitted from the laser, and the cooling systems typically accompanying laser applicators. Preliminary computations of the effect of the laser fiber active cooling system have shown substantial effects on the resulting thermal distribution.
In summary, results demonstrate that a calibration pulse prior to MRgLITT can be used to substantially increase modeling accuracy for delivery prediction. Further, results indicate the existence of a critical time interval of the calibration objective function beyond which further use of data provides diminishing returns. This observation combined with the relatively limited amount of function-gradient computations needed for convergence, Figure 9, provides an achievable strict upper bound on the computational time needed to converge to a calibrated Pennes model; this may be used as a reference for computer guided MRgLITT protocol design. However, within the perspective of a clinical setting, the overhead associated with the computation of the calibration problem needs an order of magnitude speedup. Significant algorithmic changes are needed to achieve this level of performance. One possible method, that is even more suitable for dynamic feedback, is suggested by the presented data. The smaller time windows did decrease the objective function, even though not dramatically. These relatively computationally inexpensive problems may be used as initial conditions for the full calibration problems of interest and the overhead of implementing the necessary complex data structures would be beneficial. We hope that these results help guide independent in vivo experiments and translational research on the path to realizing the model calibration aspect of computer guided LITT technology within a clinical setting.
The research in this paper was supported in part through 5T32CA119930-03 and K25CA116291 NIH funding mechanisms and the National Science Foundation under grant CNS-0540033. The authors would also like to thank the ITK , Paraview , PETSc , TAO , libMesh , and CUBIT  communities for providing truly enabling software for scientific computation and visualization. The initial 3D canine model in Figure 1 was provided by AddZero and obtained from the official Blender  model repository.
D. Fuentes, The University of Texas M.D. Anderson Cancer Center, Department of Imaging Physics, Houston TX 77030, USA.
Y. Feng, Computational Bioengineering and Nanotechnology Lab, The University of Texas at San Antonio, San Antonio, TX 78749, USA.
A. Elliott, The University of Texas M.D. Anderson Cancer Center, Department of Imaging Physics, Houston TX 77030, USA.
A. Shetty, The University of Texas M.D. Anderson Cancer Center, Department of Imaging Physics, Houston TX 77030, USA.
R. J. McNichols, BioTex, Inc., Houston TX 77054, USA.
J. T. Oden, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin TX 78712, USA.
R. J. Stafford, The University of Texas M.D. Anderson Cancer Center, Department of Imaging Physics, Houston TX 77030, USA.