Temperature-Change-Based Thermal Tomography

Abstract

Thermal properties of biological tissues play a critical role in the study of tumor angiogenesis and the design and monitoring of thermal therapies. To map thermal parameters noninvasively, we propose temperature-change-based thermal tomography (TTT) that relies on relative temperature mapping using magnetic resonance imaging (MRI). Our approach is unique in two aspects: (1) the steady-state body temperature in thermal equilibrium is not restricted to be spatially invariant, and (2) absolute temperature mapping is not required. These two features are physiologically realistic and technically convenient. Our numerical simulation indicates that a (9mm)³ tumor inside a breast phantom can be reliably depicted, assuming moderate temperature mapping accuracy of 0.5^∘C.

---

Traditional thermal tomography is based on using an infrared camera to measure the surface temperature of a tissue and then solving the inverse heat transfer problem to reconstruct the interior tissue thermal parameters [1]. Since 2002, a few papers [2–4] have discussed the possibility of using magnetic resonance imaging (MRI) to obtain the temperature distribution within the tissue and reconstruct the temperature coefficients therein. Current approaches in literature, however, have a few significant drawbacks. First of all, the theoretical model in previous work [2–4] is based on the original Pennes’ bioheat transfer equation [5], which requires absolute-temperature measurements. Yet most of the MRI-based temperature-mapping techniques, such as those based on proton resonance frequency (PRF), can only provide relative-temperature-change [6]. Second, the reconstruction algorithms in previous work [3, 4] are mathematically complicated. Last but not least, at thermal equilibrium, the in vivo temperature distribution in the region of interest (ROI) is often nonuniform. For example, the presence of a diseased tissue such as a tumor often leads to increased metabolic activity and an abnormal temperature distribution. In particular, clinical studies have suggested that the local skin temperature over a breast tumor can increase significantly (by 1–3°C) over the normal breast skin temperature [7]. This effect is often ignored in current analyses, which assume that the steady-state temperature within the ROI is uniform. To address these deficiencies, we formulate a temperature-change-based thermal tomography (TTT) method and demonstrate that TTT can potentially lead to new and noninvasive medical modalities for breast cancer imaging and other important biomedical applications.

Assuming absolute temperature within a biological tissue is described by the Pennes’ bioheat transfer equation, we have 𝜌𝐶_𝑇𝜕𝑇𝜕𝑡=∇⋅(𝜅∇𝑇)−𝜔_𝐵𝐶_𝐵(𝑇−𝑇_𝐴)+𝑄_int+𝑄_ext, (1) where 𝑇 is the absolute temperature, 𝜅 is tissue thermal conductivity, 𝑇_𝐴 is the normal body temperature, 𝜔_𝐵 is the blood perfusion rate, 𝑐_𝑇 and 𝑐_𝐵 are the specific heat of tissue and blood, respectively, and 𝑄_int and 𝑄_ext are the internal and external heat sources, respectively. In the absence of external heat sources, the absolute temperature distribution 𝑇_𝑆 at thermal equilibrium satisfies ∇⋅(𝜅∇𝑇_𝑆)−𝜔_𝐵𝐶_𝐵(𝑇_𝑆−𝑇_𝐴)+𝑄_int=0. We note that the steady-state temperature 𝑇_𝑆 may not be uniform within the ROI. However, the impacts due to nonuniform 𝑇_𝑆 can be eliminated through relative temperature measurements. In relative temperature measurements, we first obtain a baseline MRI image at thermal equilibrium, which can be accomplished in vivo and noninvasively. We then thermally excite the sample and use MRI to record the thermal processes within the biological tissue. From temperature sensitive parameters such as proton resonance frequency [6], we can extract the temperature change Δ𝑇 by comparing the MRI images taken at a later time and the baseline MRI image. Since the absolute temperature 𝑇_𝑆+Δ𝑇 satisfies (1), we can recast the Pennes’ bioheat transfer equation for the relative temperature Δ𝑇 as 𝜕Δ𝑇𝜕𝑡=∇⋅[𝛼₁∇(Δ𝑇)]−𝛼₂1Δ𝑇+𝜌𝐶_𝑇𝑄_ext, (2) where 𝛼₁ is thermal diffusivity and equals to 𝜅/𝜌𝐶_𝑇, 𝛼₂ is blood perfusion parameter and is given by 𝐶_𝐵𝜔_𝐵/𝜌𝐶_𝑇, and 𝑄_ext represents heat generated by external sources (e.g., focused ultrasound). Now we continuously perform relative temperature measurements from time 𝑡₀ and to time 𝑡₁. Since Δ𝑇 satisfies (2), after integrating both sides of (2) from 𝑡₀ to 𝑡₁ in the absence of external sources, we obtain ∇⋅(𝛼₁∇⟨Δ𝑇⟩)−𝛼₂⟨Δ𝑇⟩=𝑇(𝑡₁)−𝑇(𝑡₀)Δ𝑡, (3) where Δ𝑡=𝑡₁−𝑡₀, and ⟨Δ𝑇⟩ is the time average of the relative temperature, that is, ⟨Δ𝑇(𝑥,𝑦,𝑧)⟩≡∫^𝑡1_𝑡0Δ𝑇(𝑥,𝑦,𝑧,𝑡)𝑑𝑡/Δ𝑡. If we discretize the spatial domain of the biological structure into 𝑁 grid points and apply (3) at each point, we obtain 𝑁 linear equations with 2𝑁 unknown variables (one 𝛼₁ and one 𝛼₂ at each grid point). Performing an additional set of MRI temperature measurements with a different initial temperature, we have 2𝑁 linear equations and can uniquely determine the values of 𝛼₁ and 𝛼₂ over the entire numerical grid. Such a direct approach for TTT is significantly simpler than other thermal reconstruction methods in literature [2–4].

Having established the general framework of TTT, we consider an important practical application: the imaging and detection of breast cancer. It has been shown that blood flow in tumors and in normal tissues can differ dramatically [1, 2]. For example, the blood perfusion rate ?_B of normal breast tissue is 0.7?kg⋅m⁻³⋅s⁻¹ [8], whereas 𝜔_𝐵 for breast adenocarcinoma is 5∼6.7kg⋅m⁻³⋅s⁻¹ [2]. Assuming a tissue density 𝜌 of 1×10³kg⋅m⁻³, a tissue specific heat of 𝑐_𝑇=3550J⋅kg⁻¹⋅^∘C⁻¹ [1, 9], and a blood specific heat of 𝑐_𝐵=3600J⋅kg⁻¹⋅^∘C⁻¹ [1], the blood perfusion parameter 𝛼₂ for normal breast is 7.1×10⁻⁴s⁻¹, whereas for breast tumor, 𝛼₂ ranges from 5.1×10⁻³s⁻¹ to 6.8×10⁻³s⁻¹. Notice that the value of 𝛼₂ in cancerous tissue increases by almost one order of magnitude. Thermal conductivity 𝜅, however, remains similar in cancerous and normal tissues. For example, the values of 𝜅 for tumor, lung, and skeletal muscle are 0.642 [9], 0.30∼0.55 [1], and 0.45∼0.55W⋅m⁻¹⋅^∘C⁻¹ [1], respectively. In this letter, we assume that 𝜅 is 0.47?W⋅m⁻¹⋅^∘C⁻¹ for both the tumor and the normal tissue, which gives a thermal diffusivity 𝛼₁ of 0.13?mm²⋅s⁻¹. With these considerations, a breast phantom is established in Figure 1: the phantom is a hemisphere with a radius of 75?mm and contains a (9?mm)³ tumor buried within. The thermal diffusivity 𝛼₁ is 0.13?mm²⋅s⁻¹ throughout the phantom (including the tumor). The blood perfusion parameter 𝛼₂ is 5.6×10⁻³s⁻¹ within the tumor and 7.1×10⁻⁴s⁻¹ outside of it.

Figure 1 

Schematic of the breast phantom.

Now we use a finite-difference (FD) method [10] to simulate temperature evolution within the breast phantom. (Later on, the simulated temperature field is used for thermal reconstruction tests). The simulation domain, which is illustrated in Figure 1, has a dimension of 150×150×84?mm³. It is evenly divided into an array of FD cells with a cell size of Δ𝑥=Δ𝑦=Δ𝑧=3?mm. The number of cells along the 𝑥, 𝑦, and 𝑧 direction is, respectively, 50, 50, and 28. Each FD cell is labeled by three indices as (𝑖, 𝑗, 𝑘), with cell center located at the position of (𝑖Δ𝑥, 𝑗Δ𝑦, 𝑘Δ𝑧). The center of the breast phantom is at cell (25, 25, 28). The tumor region is represented by 9 FD cells whose indices are in the range of 27≤𝑖≤29, 19≤𝑗≤21, and 18≤𝑘≤20. The FD simulation is carried out by discretizing (2) and evolving the temperature field forward with a time step of Δ𝑡=1 second. For boundary condition, we set the relative temperature shift at the breast boundary to be zero. This assumption is reasonable, since the temperature at the breast surface can be kept at a constant value and deep within the breast (at the top 𝑧 surface) the temperature shift should be minimal. For the initial temperature at time 𝑡₀, we assume that it has a Gaussian profile Δ𝑇(𝑥,𝑦,𝑧)=Δ𝑇₀exp{−[(𝑥−𝑥₀)²+(𝑦−𝑦₀)²+(𝑧−𝑧₀)²]/𝛿²} with the following parameters: Δ𝑇₀=5^∘C, 𝑥₀=25Δ𝑥, 𝑦₀=24Δ𝑦, 𝑧₀=17Δ𝑧, and 𝛿=9Δ𝑥. In Figure 2, we show a few snapshots of the relative temperature field obtained through FD simulations. The results in Figure 2(a) are plotted along a straight line that passes through the tumor center and is along the 𝑥 direction, that is, the dashed line in Figure 1. Since MRI has a finite accuracy in temperature measurement, to make our analysis more realistic, we add Gaussian noises with a standard deviation of 0.5°C to the simulation results. (A measurement accuracy of 0.5°C with a spatial resolution of ∼1?mm is realistic given current MRI technology [2]). Figure 2(b) gives the temperature field at 𝑡=200?s, one from FD simulation and the other with noise added. Figure 2(c) shows the temporal evolution of the noise-added temperature field at cell (28, 20, 19), the center of the tumor.

(a)

(b)

(c)

Figure 2 

(a) Relative temperature field plotted along the dashed line in Figure 1, at time 𝑡=0, 𝑡=200 seconds, 𝑡=400 seconds, respectively. (b) Relative temperature field at 𝑡=200 seconds. The dots and the crosses represent FD simulation data and noise-added data, respectively. (c) Time evolution of the noise-added temperature field at FD cell (28, 20, 19). The dashed line is an exponential fit.

We now treat the noise-added simulation results as MRI measurements data and use them to test our thermal reconstruction algorithm. In this letter, we assume that 𝛼₁ is known to be 0.13?mm²⋅s⁻¹ throughout the breast phantom, which means that we only need one set of temperature measurements to reconstruct the blood perfusion parameter 𝛼₂. This assumption is reasonable, since the value of 𝛼₁ does not vary significantly from tumor to normal tissues. Furthermore, the importance of the thermal diffusion term can be greatly reduced by choosing a relatively flat temperature distribution, which is the case for the initial Gaussian field used in Figure 2. To further illustrate this point, consider an extreme case: for a uniform temperature distribution, the thermal diffusion term in (2) becomes zero, and the value of 𝛼₁ has no impact on the temperature decay rate.

In the absence of noise, the spatial distribution of 𝛼₂ can be directly obtained through (3). However, if the input temperature-change data Δ𝑇 contain significant noises, a straightforward application (3) is not advisable. In fact, a finite-difference evaluation of ∇⋅(𝛼₁∇⟨Δ𝑇⟩) can significantly amplify the already nonnegligible measurement noises. Fortunately, we can use numerical fitting to achieve significant noise reduction. As an example, we notice that with a uniform distribution of 𝛼₁, ∇⋅(𝛼₁∇⟨Δ𝑇⟩) becomes 𝛼₁∇²⟨Δ𝑇⟩, which can be evaluated by polynomial fitting. More specifically, to evaluate 𝜕²⟨Δ𝑇⟩/𝜕𝑥² at cell (𝑖₀,𝑗₀,𝑘₀), we use the polynomial ⟨Δ𝑇⟩=𝑎₂(𝑖−𝑖₀)²Δ𝑥²+𝑎₁(𝑖−𝑖₀)Δ𝑥+𝑎₀ to fit the temperature data at cells (𝑖,𝑗₀,𝑘₀), with 𝑖 in the range of 𝑖₀−2≤𝑖≤𝑖₀+2. After fitting, the value of 𝜕²⟨Δ𝑇⟩/𝜕𝑥² is simply given by the fitting coefficient 𝑎₂. 𝜕²⟨Δ𝑇⟩/𝜕𝑦² and 𝜕²⟨Δ𝑇⟩/𝜕𝑧² can be similarly evaluated. Next, we consider the evaluation of 𝑇(𝑡₁)−𝑇(𝑡₀) in (3). To eliminate measurement noises at any specific time, we use an exponential function to fit the temporal evolution of the temperature field. In Figure 2(c), we show the fitting result (represented by a dashed line) at cell (28, 20, 19). (The original noise-added temperature data are represented as dots.) As can be seen from Figure 2(c), the exponential fitting provides a reasonable approximation to the noise-added temperature data. We can then use the fitted curve to calculate 𝑇(𝑡₁)−𝑇(𝑡₀). After applying the aforementioned fitting schemes to (3), we can reconstruct the spatial distribution of the blood perfusion parameter 𝛼₂ from the noise-added FD simulation data. The reconstruction results are shown in Figure 3 as diamonds. (The data are plotted along the dashed line in Figure 1, which is the same as the line used in Figure 2). We repeat the same reconstruction procedure 30 times, each time using the same FD simulation data with randomly generated Gaussian noises. From the statistical fluctuations of the reconstructed 𝛼₂, we can obtain the standard deviations for 𝛼₂, which are given in Figure 3 as error-bars. We also apply the same reconstruction procedure to the FD simulation data without any noise. The results are shown in Figure 3 as dots connected by solid lines. The input 𝛼₂ values (used in forward FD simulations) are also shown in Figure 3 as a dashed line. From Figure 3, we find that despite the 0.5°C random noises, we can successfully reconstruct the spatial distribution of the blood perfusion parameter 𝛼₂. We also observe that the numerical fitting introduces approximately 20% systemic errors to the thermal reconstruction. However, it is possible to reduce the systemic error by developing a more accurate noise reduction algorithm.

Figure 3 

Reconstructed thermal coefficients 𝛼2 along the dashed line in Figure 1.

We emphasize that the advantages of TTT are derived primarily from its relative-temperature basis. By formulating a differential equation based on temperature changes (i.e., (2), we remove any potential impact due to an abnormal temperature profile in the ROI at thermal equilibrium. Furthermore, the requirement of relative temperature measurements makes TTT much more suitable for clinical applications.

The algorithm presented in this paper serves as a proof-of-concept validation and can be improved in future studies. For example, the method of five-point polynomial fitting can be replaced by other spatial filtering techniques that can effectively reduce the measurement noises. The exponential temporal fitting can also be replaced by more sophisticated filtering techniques such as the Wiener filter. The assumption of an initial Gaussian temperature distribution is also not necessary: We have verified the reconstruction algorithm (i.e., (3) using different initial temperature distribution and other means of thermal excitations.

Compared with other methods of thermal tomography in current literature, the method presented in this paper is significantly simpler. For example, the method presented in [3] relies on the framework of Bayesian statistics and statistical inverse problems. Yet in current work, we can achieve similar resolution using a much simpler numerical approach. We also point out that within a temperature range of less than 4°C (refer to Figure 2), tissue thermal parameters should remain approximately the same [9]. Consequently, the effects of nonlinear heat transfer should not significantly change the main conclusions of this work. As presented in current paper, TTT may require an image acquisition time of a few minutes, during which the effects of motion artifacts may arise. Such artifacts, however, can be addressed using methods reviewed in [6].

In conclusion, we have developed a temperature-change-based thermal tomography and established a simple and straightforward thermal reconstruction algorithm. Using realistic thermal parameters and assuming a temperature measurement accuracy of 0.5°C, we have successfully reconstructed the spatial distribution of a (9?mm)³ tumor embedded in a human breast phantom. We believe that TTT can lead to new and noninvasive medical modalities for cancer detection and imaging.

Acknowledgments

The authors gratefully acknowledge the support of the Institute for Critical Technology and Applied Science (ICTAS) at Virginia Tech. They would also like to thank Professor Robert Kraft of Wake Forest University Health Sciences for helpful discussions on MRI.