Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3857613

Formats

Article sections

- Abstract
- I. Introduction
- II. The Calibration Problem
- III. Experimental Setup
- IV. Results
- V. Discussion and Conclusion
- References

Authors

Related links

IEEE Trans Biomed Eng. Author manuscript; available in PMC 2013 December 10.

Published in final edited form as:

Published online 2010 February 5. doi: 10.1109/TBME.2009.2037733

PMCID: PMC3857613

NIHMSID: NIHMS526379

D. Fuentes, The University of Texas M.D. Anderson Cancer Center, Department of Imaging Physics, Houston TX 77030, USA;

D. Fuentes: dtfuentes/at/mdanderson.org; Y. Feng: yusheng.feng/at/utsa.edu; A. Elliott: andrew.elliott/at/mdanderson.org; A. Shetty: anil.shetty/at/mdanderson.org, ; R. J. McNichols: roger/at/biotexmedical.com; J. T. Oden: oden/at/ices.utexas.edu

The publisher's final edited version of this article is available at IEEE Trans Biomed Eng

See other articles in PMC that cite the published article.

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 [9]. The cyberinfrastructure [1], [6] 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 [4] and substantial translational research and validation is needed to fully realize the potential of this technology [18], [21] 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 [7]. 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.

Schematic diagram of *in vivo* MR-guided LITT calibration procedure in a canine model of prostate. Contrast enhanced T1-W MR images have been volume rendered to better visualize the relationship of the target volume and applicator trajectory to the surrounding **...**

The problem of bioheat transfer model calibration is to determine the set of thermal parameters that minimize the *L*_{2}(0, *T; L*_{2}(Ω)) norm of the difference between the predicted temperature field, *u*(**x**, *t*) and the temperature field observed *in vivo* thermal images of the experiment, *u ^{MRTI}*(

(1)

*dx* = *dx*_{1}*dx*_{2}*dx*_{3} 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 [16] of bioheat transfer with a isotropic laser heat source.

Pennes model has been shown to provide very accurate prediction of bioheat transfer [5], [14], [23] and is used as the basis of the finite element prediction. The full initial boundary value model is defined by the following system:

(2)

The initial temperature field, *u*(**x**, 0) = *u*^{0}, 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}*,

where
, *k*_{3} [*K*] . Perfusion is modeled with a nonlinear dependence on temperature, Figure 3:

where
, *ω _{NI}* [

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 [3] 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 [15]

Here *p* is the solution to the adjoint problem and _{0}(**x**), _{1}, _{2}, _{3}, _{0}(**x**), * _{N}*,

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

(a) A cutplane though the FE prostate model illustrating the post-calibration temperature prediction of Pennes model. The cutplane shown is intersecting and coplanar with the catheter. The temperature is given in °C. Mesh lines illustrate the **...**

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) [11] 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 [19]. 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 [11]) 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

Time varying plots of the post-calibration pointwise error in the model prediction, *u* – *u*^{MRTI}, along the x-axis of Figure 4(a) and of the thermal image data along the y-axis of Figure 4(a). The temperature and distance are given in units of Celsius **...**

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 comparison of cutplanes through the (a) measured MR temperature image and (b) the pre and post calibration finite element prediction is shown. The cutplanes are taken at the same axial position within the anatomy; the plane is 3mm from the plane that **...**

The top, middle, and bottom cutlines Figure 6(a) are used to compare the pointwise thermal image temperature and pre and post calibration finite element prediction of the temperature versus cutline distance. The profiles illustrate the effect of the heterogeneity. **...**

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.

A summary of a calibration study using various time windows of thermal image information is shown. The power versus time profile of the calibration pulse is shown in the insert; the pulse is held constant at 5W for Δ*t* =60 seconds then turned off. **...**

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-

The number of function evaluations required for convergence of the Quasi-Newton optimization method used in the calibration study is shown. The graph is intended to convey the general convergence behavior observed from the breadth of the calibration problems **...**

(3)

For a particular problem, when *Q _{converged}* 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 [9], 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 [12], Paraview [10], PETSc [2], TAO [3], libMesh [13], and CUBIT [20] 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 [17] model repository.

Webpage: http://wiki.ices.utexas.edu/dddas

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.

1. Atkins DE, Droegemeier KK, Feldman SI, Garcia-Molina H, Klein ML, Messerschmitt DG, Messina P, Ostriker JP, Wright MH. Technical report. NSF; 2003. Revolutionizing science and engineering through cyberinfrastructure: Report of the national science foundation blue-ribbon advisory panel on cyberinfrastructure.

2. Balay Satish, Gropp William D, McInnes Lois C, Smith Barry F. Technical Report ANL-95/11 - Revision 2.1.5. Argonne National Laboratory; 2003. Petsc users manual.

3. Benson Steven J, McInnes Lois Curfman, Moré Jorge, Sarich Jason. Technical Report ANL/MCS-TM-242. Mathematics and Computer Science Division, Argonne National Laboratory; 2005. TAO user manual (revision 1.8) http://www.mcs.anl.gov/tao.

4. Carpentier A, McNichols RJ, Stafford RJ, Itzcovitz J, Guichard JP, Reizine D, Delaloge S, Vicaut E, Payen D, Gowda A, et al. Real-time magnetic resonance-guided laser thermal therapy for focal metastatic brain tumors. Neurosurgery. 2008;63(1 Suppl 1):8. [PubMed]

5. Charny CK. Mathematical models of bioheat transfer. Adv Heat Trans. 1992;22:19–155.

6. Darema F. DDDAS: Dynamic Data Driven Applications Systems. website: www.nsf.gov/cise/cns/dddas.

7. 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. Numerical Implementation of Bioheat Models and Equations, chapter 9: Computational Infrastructure for the Real-Time Patient-Specific Treatment of Cancer. Vol. 3. Taylor & Francis Group; 2008. Advances in Numerical Heat Transfer.

8. Diller KR, Valvano JW, Pearce JA. Bioheat transfer. In: Kreith F, Goswami Y, editors. The CRC Handbook of Mech Engr. 2. CRC Press; Boca Raton: 2005. pp. 4–278–4–357.

9. Fuentes D, Oden JT, Diller KR, Hazle J, Elliott A, Shetty A, Stafford RJ. Computational modeling and real-time control of patient-specific laser treatment cancer. Ann BME. 2009;37(4):763. [PubMed]

10. Henderson A, Ahrens J. The ParaView Guide. Kitware; 2004.

11. Hindman JC. Proton Resonance Shift of Water in the Gas and Liquid States. The Journal of Chemical Physics. 1966;44:4582–4592.

12. Ibanez L, Schroeder W, Ng L, Cates J. The ITK Software Guide. 2. Kitware, Inc; 2005. http://www.itk.org/ItkSoftwareGuide.pdf.

13. Kirk BS, Peterson JW. libMesh-a C++ Finite Element Library. CFDLab; 2003. URL http://libmesh.sourceforge.net.

14. Liu J, Zhu L, Xu L. Studies on the three-dimensional temperature transients in the canine prostate during transurethral microwave thermal therapy. J Biomech Engr. 2000;122:372–378. [PubMed]

15. 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. Num Meth PDE. 2007;23(4):904–922. [PMC free article] [PubMed]

16. Pennes HH. Analysis of tissue and arterial blood temperatures in the resting forearm. J Appl Physiol. 1948;1:93–122. [PubMed]

17. Roosendaal T, Wartmann C. The official Blender 2.0 guide. Prima Tech; 2000.

18. Salomir R, et al. Hyperthermia by MR-guided focused ultrasound: accurate temperature control based on fast MRI and a physical model of local energy deposition and heat conduction. Magn Reson Med. 2000;43(3):342–347. [PubMed]

19. Stafford RJ, Price RE, Diederich CJ, Kangasniemi M, Olsson LE, Hazle JD. Interleaved echo-planar imaging for fast multiplanar magnetic resonance temperature imaging of ultrasound thermal ablation therapy. Journal of Magnetic Resonance Imaging. 2004;20(4):706–714. [PubMed]

20. Blacker T, et al. Cubit Users Manual. 2008 http://cubit.sandia.gov/documentation.

21. Vimeux FC, et al. Real-time control of focused ultrasound heating based on rapid MR thermometry. Invest Radiol. 1999;34(3):190–193. [PubMed]

22. Welch AJ, van Gemert MJC. Optical-Thermal Response of Laser-Irradiated Tissue. New York: Plenum Press; 1995.

23. Wissler EH. Pennes’ 1948 paper revisited. J Appl Physiol. 1998;85:35–41. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |