Non-invasive tissue temperature measurements based on quantitative diffuse optical spectroscopy (DOS) of water

We describe the development of a non-invasive method for quantitative tissue temperature measurements using Broadband diffuse optical spectroscopy (DOS). Our approach is based on well-characterized opposing shifts in near-infrared (NIR) water absorption spectra that appear with temperature and macromolecular binding state. Unlike conventional reflectance methods, DOS is used to generate scattering-corrected tissue water absorption spectra. This allows us to separate the macromolecular bound water contribution from the thermally induced spectral shift using the temperature isosbestic point at 996 nm. The method was validated in intralipid tissue phantoms by correlating DOS with thermistor measurements (R = 0.96) with a difference of 1.1 ± 0.91 °C over a range of 28–48 °C. Once validated, thermal and hemodynamic (i.e. oxy- and deoxy-hemoglobin concentration) changes were measured simultaneously and continuously in human subjects (forearm) during mild cold stress. DOS-measured arm temperatures were consistent with previously reported invasive deep tissue temperature studies. These results suggest that DOS can be used for non-invasive, co-registered measurements of absolute temperature and hemoglobin parameters in thick tissues, a potentially important approach for optimizing thermal diagnostics and therapeutics.


Introduction
Temperature-dependent changes in near-infrared (NIR) water absorption spectra have been well characterized (Collins 1925, Otal et al 2003. As temperature increases, the water peak around 970 nm blue-shifts, narrows and increases in intensity (Collins 1925, Kelly et al 1995, Libnau et al 1994, McCabe et al 1970, Merritt 2005, providing a sensitive spectral signature of temperature. Since tissue is primarily composed of water (typically ∼60-70%), tissue water spectra can potentially reveal tissue temperature. However, at depths greater than the NIR transport scattering length in tissues (e.g. 1 mm), multiple light scattering obscures temperature-dependent water spectral shifts, impairing reliable spectral measurements of tissue temperature.
Broadband diffuse optical spectroscopy (DOS) acquires scattering-separated absorption spectra, typically from 650 to 1000 nm. In this work, we employ a type of DOS that combines broadband frequency-domain photon migration (FDPM) with steady-state (SS) reflectance spectroscopy (Tromberg et al 1993, Bevilacqua et al 2000, Cerussi et al 2006. This approach provides high resolution, quantitative absorption spectra, which are needed to characterize temperature-dependent changes associated with the water absorption features. DOS absorption spectra are used to quantify major tissue NIR chromophores including oxy-and deoxy-hemoglobin, lipid and water at centimeter depths. By measuring scattering-corrected absorption changes in the 970 nm water peak, additional information regarding the disposition of tissue water can be obtained. This includes the bound water fraction (Chung et al 2008) and, as we demonstrate in this paper, tissue temperature.
DOS measurements of temperature, hemodynamic and metabolic parameters are colocalized. As a result, this non-invasive technique could potentially be employed to optimize dosimetry during thermotherapies, enhancing therapeutic outcome and reducing damage to normal tissues. Thermal therapies are widely used in medicine, spanning from cryosurgery to hyperthermia and thermal ablation. Non-invasive imaging methods that measure and map subsurface tissue temperature can be difficult and expensive to employ (e.g. MRS), and implantable sensors may perturb energy deposition and tissue temperature distribution. Methods such as hyperthermia utilize heat to remove tumors by killing cells directly or sensitizing them to chemotherapy agents and/or radiation (Kim and Hahn 1979, Rieke and Pauly 2008, Seegenschmiedt et al 1995, Thomsen 1991, Vargas et al 2004. In order to remove a tumor using hyperthermia, heat in the range of 43-50 • C is applied for an extended period, in some cases hours (Vargas et al 2004). Vargas et al showed that tumor necrosis fraction is determined by peak tumor temperature as well as the thermal dose, and a 2.6 • C difference in the peak tumor temperature can result in a 50% difference in the tumor necrosis. Thus, absolute temperature measurements could be used to optimize dosimetry while co-localized tissue functional parameters such as oxy-and deoxy-hemoglobin report on the biologic effects of therapy.
In this work we validate our approach for quantifying temperature in multiple-scattering intralipid phantoms and compare DOS temperature calculations to simultaneous thermistor measurements. Our results show a correlation between non-invasive optical measurements and invasive thermal measurements (R = 0.96, average difference = 1.1 ± 0.91 • C, in 28-48 • C). In addition, we demonstrate the potential for broadband DOS to measure absolute temperature and hemodynamic parameters in vivo in human subjects during mild cold stress. Pure water absorption spectra measured as a function of temperature (15-65 • C). The water peak around 970 nm increases in intensity, narrows in width and blue-shifts in wavelength as temperature increases. An isosbestic point at 996 nm is observed.

Background theory
A pronounced water absorption peak appears in the NIR around 970 nm. This peak is due to the combination of the first harmonic of the O-H symmetric stretch vibration and the fundamental anti-symmetric stretch vibration from hydrogen bound O-H (Bayly et al 1963).
The characteristics of this peak are sensitive indicators of the local environment of water molecules. The spectral features of water absorption depend upon the temperature and water binding state. As temperature increases, the fraction of hydrogen-bound water molecules is reduced, causing the 970 nm water peak to increase in intensity, narrow in bandwidth and shift to higher energy (i.e. blue-shift), as shown in figure 1 (Kelly et al 1995, Libnau et al 1994, McCabe et al 1970, Merritt 2005. Walrafen (1967) explained that these changes are from a shift in the equilibrium between hydrogen-bound and -unbound (free) water molecules. Hydrogen bonding between water and macromolecules such as proteins causes further spectral broadening and red shifting of the 970 nm water peak (Pimentel andMcClellan 1960, Chung et al 2008). These two competing effects (i.e. red and blue-shifts of the peak) complicate the proper interpretation of the water local environment. However, there is also an isosbestic point with respect to temperature in the water absorption spectrum at 996 nm (figure 1, Buijs and Choppin 1963, Walrafen 1967, Hollis 2002, Merritt 2005. As a result, the absorption coefficient at 996 nm can be used to account for macromolecular bound water effects independent of tissue temperature. Once the bound water correction is determined, the best fit to the temperature-dependent pure water spectra yields the absolute temperature of the measured tissue volume.

Instrumentation
Broadband DOS employed in this study combines multi-frequency, frequency domain photon migration (FDPM) with broadband steady-state (SS) spectroscopy in order to acquire continuous NIR absorption and scattering spectra in turbid media (Jakubowski et al 2009). Multi-frequency FDPM separates the effects of absorption from scattering in tissues using model-based analysis Gratton 1993, Tromberg et al 1993). Models include a P1 approximation to the transport equation in the semi-infinite (for in vivo measurements) or infinite (for phantoms) geometry using an extrapolated boundary condition (Fishkin et al 1996, Haskell et al 1994. The SS component enables the acquisition of continuous absorption spectra of more than 1000 wavelengths from 650 to 998 nm for phantoms and 1013 nm for in vivo tissues (B&W Tek, Newark, DE, model 611a and 611b, respectively). Because the spectrometer is auto-calibrated, the absolute wavelengths are obtained for every measurement. Instrument response is determined using tissue-simulating phantoms with known scattering and absorption values for FDPM and spectraflect-coated integrating spheres with known reflectance values for SS measurements, as described in Cerussi et al (2006). Details of the broadband DOS system have been previously reported (Bevilacqua et al 2000, Cerussi et al 2006, Jakubowski et al 2009. Measurements were performed by placing fixed source and detector optical fibers in the middle of intralipid phantoms (infinite geometry) or on a forearm (semi-infinite geometry) of a human subject. The multi-frequency (∼400 frequencies from 50 to 300 MHz) swept FDPM data were acquired for 661, 681, 783, 806, 823 and 850 nm lasers (<20 mW power) for the phantoms and 658, 682, 785, 810, 830 and 850 nm lasers for the in vivo tissue measurements. The SS data (single tungsten-halogen lamp) for all wavelengths with 0.3-0.5 nm increments between 650 and 998 (phantoms) or 1013 nm (in vivo) were acquired in a total of ∼6 s. The differences between the two systems were due to the availability of components at different times. The source-detector separation was 9 mm for phantom measurements and 20 mm for in vivo measurements. Model fits to SS and FDPM data were performed in order to recover absorption (μ a ) and reduced scattering (μ s ) values for all wavelengths between 650 and 998 nm as described in Bevilacqua et al (2000), Cerussi et al (2006), and Jakubowski et al (2009).

Spectral processing algorithm
The extinction coefficient spectra of major chromophores in tissues and phantoms (oxy-and deoxy-hemoglobin, water and lipid for tissues; water, lipid and nigrosin in phantoms) were used to recover the concentrations of these components from measured absorption spectra (Cerussi et al 2002, 2007, Cubeddu et al 2000, Pogue et al 2004, Quaresima et al 1998. For pure water, distilled water placed in a cuvette was measured using a spectrophotometer (Beckman DU 650) at various temperatures (Merritt 2005). Nigrosin (Sigma-Aldrich, St. Louis, MO, water soluble) was also measured with the spectrophotometer using a 511 mg l −1 stock solution. The lipid spectra were obtained from mammalian fat by van Veen et al (2005). For oxy-and deoxy-hemoglobin, we employed molar extinction coefficients of Zijlistra et al (2000). The measured absorption spectrum and extinction coefficient spectra provided needed information to solve nonnegative least-squares constraints problems (using 'lsqnonneg' in MATLAB R ) to recover the concentrations of the major absorbers in phantoms and tissues mentioned above. A 'water only' tissue absorption spectrum was obtained by eliminating contributions from all non-water tissue absorbers. This was accomplished by subtracting hemoglobin (oxy-and deoxy-) and lipid extinction coefficient spectra, multiplied by their concentrations, from broadband tissue absorption spectra (Chung et al 2008). The 'water only' tissue absorption spectrum was normalized by multiplying the ratio of the water peak of a pure water spectrum to the measured water peak. The pure water spectrum used at this step was of the temperature that was measured using a thermistor for the phantoms or 33.5 • C for the forearm measurement.
In order to estimate temperature using the water absorption spectrum, we employed a twostep process. First, we corrected the tissue water spectra for the spectral shift and broadening caused by macromolecular binding using an iterative fitting algorithm and the isosbestic point of the pure water spectra. Then, the best-fit pure water spectrum at a specific temperature to the bound water corrected spectrum was used to calculate tissue temperature. These steps are explained in detail below.
Step 1. Bound effect correction In order to determine the spectral shift in the absorption spectrum due to bound water, we fit the extinction coefficient spectrum for each tissue component to the absorption spectrum while correcting the water extinction coefficient spectrum using an iteration algorithm. A pure water extinction coefficient spectrum at the pre-measured (for phantoms) or assumed (for tissues) temperature is shifted using equation (1) (Merritt 2005): where wpwl is the peak wavelength of the pure water spectrum, bws is the bound water shift defined as the peak-to-peak distance between the tissue water spectrum and the corrected pure water spectrum, and WL is the upper wavelength limit used for the bound water shift correction. It was found that above 935 nm the tissue water absorption spectrum was linearly shifted to higher wavelengths relative to the pure water spectrum; thus, 935 nm was used as the starting wavelength for the bound water shift calculation. The chromophore fit is then executed with the newly corrected pure water spectrum and a loop is entered in which the chromophore fit and the bound water shift calculation are performed until the bws in equation (1) changes by less than 0.1 nm. The calculated bws in equation (1) was used in order to blue-shift a normalized tissue water spectrum as the first step of bound water effect correction. Although significant bound water shifting is accounted for using the above step, absorption values at the isosbestic point (996 nm, shown in figure 1) are not fully corrected. This may be due to uncertainties in our estimation of the temperature of pure water basis spectra with respect to tissue water. Additional correction is therefore performed on the normalized spectra by forcing an isosbestic point at 996 nm. We validated this procedure by determining the relationship between the normalized intensity at 996 nm and the correction amount using two intralipid phantoms with known temperature (28-48 • C with 3 ± 1.3 • C increment) and observed a linear relationship (y = 0.7238x − 0.0303, R 2 = 0.9681) (figure 2). This equation was then used for the final blue-shifting of the tissue water spectrum in order to minimize the difference between measured and ideal spectra from 925 nm to 996 nm, where 925 nm is essentially independent of bound water effects. During the blue-shifting, the spectrum was extrapolated slightly if the upper limit of the spectrum was less than 998 nm.
Step 2. Spectral fitting for absolute temperature determination After the bound water effect corrections, the remaining tissue water spectrum was fit to a library of pure water spectra at different temperatures using a nonlinear least-squares curvefitting algorithm ('lsqcurvefit' in MATLAB R ). The temperature of the pure water spectrum that fits the tissue water spectrum with the lowest residuals was selected as the temperature of the measured tissue volume. For this fitting process, a spectral window from 940 to 995 nm was used because empirically we found it to be the most robust range for temperature fitting.

Phantom and human subjects measurements
Intralipid tissue phantoms were used for comparing broadband DOS and thermistor measurements of temperature. 600 ml of distilled water was mixed with 150 ml Intralipid (Intralipid R 20% I.V. emulsion, Fresenius Kabi, Uppsala, Sweden, distributed by Baxter, IL) and 5 ml nigrosin (Sigma-Aldrich, St Louis, MO) of 511 mg l −1 stock solution. Four of the intralipid phantoms were made and measured on different dates. A calibrated thermistor with ±0.2 • C accuracy (Newport, Santa Ana, CA) measured the temperature of the phantoms. The temperature in the middle of the phantom was measured simultaneously by the broadband DOS instrument and the thermistor. So as not to interfere with the optical measurement, the thermistor was placed about 5 mm longitudinally lower and 5 mm transversely away from the end of the detector fiber and the optical field of view. Phantom temperature was controlled using a hot plate placed underneath the phantom.
For in vivo tissue measurements, forearms (brachioradial muscles) of four healthy male subjects (20.5 ± 1.3 years) were measured continuously over a 70 min period. Ultrasound was used to measure adipose layer thicknesses, which ranged from 2 to 4 mm. The patient was lying in a supine position during the measurements. The arm was relaxed and put in a water tank that was lower than the bed height. The DOS probe was fixed onto the arm in order to maintain consistent contact with the tissue. The temperature of water in the tank was 37 ± 1 • C for the initial 5 min, decreased to 15 ± 1 • C for 30 min (cooling phase) and then increased to 37 ± 1 • C for 30 min (recovery phase) for three subjects. One of the subjects had a shorter cooling phase and recovery phase (20 min each). The temperature of water in the tank was monitored with a thermometer and kept constant by adding warm water or ice for the 37 ± 1 • C and 15 ± 1 • C water, respectively.

Results
Representative bound water and temperature-corrected water spectra from a human subject (forearm muscle) are shown in figure 3. As a result of the first bound water shift correction, Bound water corrected tissue water spectrum fit to a pure water spectrum at a specified temperature. The temperature of the best-fit pure water spectrum is the temperature of the measured tissue volume. the water absorption spectrum was blue-shifted ('first bw corrected' in figure 3(a)). The calculated bws used for this blue-shift in Step 1 was 3-6 nm for the forearm measurements depending on the subject. The bws for the intralipid phantoms was 0 nm; thus, Step 1 was skipped. The value used for the blue-shift in Step 1 for a set of longitudinal measurements is the minimum bws obtained during the baseline or initial measurements. The 'second bw corrected' spectrum in figure 3(a) is the result of the additional optimization of the blue-shift. The additional shift resulted in the best fit to one of the temperature-dependent pure water spectra as shown in figure 3(b). The temperature of the best-fit pure water spectrum was assumed to be the average temperature of the measured volume.
Temperature changes in the intralipid phantoms were monitored using both broadband DOS and a thermistor (figure 4). Temperature was measured in 0.27 ± 0.25 • C increments  during the heating process in four phantoms and during the cooling process in two of them. No significant differences between heating and cooling measurements were observed. A high correlation between the thermistor and broadband DOS measured temperature was obtained (R = 0.96). The slope of the linear fit to the temperature from both techniques was 0.94. The overall mean difference between DOS and thermistor measurements was 1.1 ± 0.91 • C. The standard error of the DOS-measured temperature was 1.0 ± 0.74 • C. Thermal changes in the human forearm were measured longitudinally in four male subjects (age: 20.5 ± 1.3). In figure 5, one of the subject's thermal and hemodynamic responses to environmental temperature changes are shown. During the 5 min baseline measurement, the temperature did not change significantly (27.6 ± 0.6 • C). Shortly after the arm was immersed in 15 • C water, arm temperature slowly decreased to 19.5 • C, about 8 • C lower than the baseline. Fifteen minutes after replacing the cold water with 37 • C water, the arm temperature recovered nearly to the initial values. In this subject, oxyhemoglobin, a sensitive indicator of tissue oxygen delivery, dropped significantly during cold stress. Deoxyhemoglobin, which reflects tissue oxygen consumption, did not change significantly. In all of the subjects, the baseline temperature during the first 5 min and the recovered temperature after the cold stress (last 5 min) were 27.6 ± 0.5 • C and 27.2 ± 0.3 • C, respectively (table 1). The minimum temperature during the cold stress was 22.8 ± 2.6 • C. The dynamic range of thermal response during the recovery was similar to the thermal dynamic range during the cold stress (0.63 ± 0.3 difference). The hemodynamic rate was different in all subjects in response to the environmental temperature changes.

Discussion
We have described a method for absolute temperature measurements based on resolving water vibrational frequency shifts due to macro-molecular binding. Because this approach accounts for thermally induced changes in light scattering, protein denaturation/aggregation and/or changes in lipids and molecular complexes that typically occur with heating do not impact accuracy. The spectral changes are observed as shifts and broadening/narrowing of the water absorption peak at 970 nm as shown in figure 1. Because the method requires quantifying subtle changes in the NIR water absorption features, efficient separation of absorption from scattering and a high-resolution calibrated spectrometer are required. The spectral processing algorithm separates counteracting bound water effects from temperature effects on the NIR water absorption spectrum. The bound water correction is based on our previous work (Chung et al 2008) where a bound water index (BWI) was introduced. The residual between a tissue water absorption spectrum and a pure water absorption spectrum was used to calculate the BWI, an index that represents both the water peak shift and its shape change. Those features were used to correct for the bound water effect on the water spectrum in this study. After the first correction using the fitting and iteration algorithm, the second correction using the pure water temperature isosbestic point at 996 nm completed the bound water effect correction (figures 2 and 3(a)). The absorption at 925 nm of normalized water absorption spectra was 0.0108 ± 0.0007 mm −1 in the intralipid phantoms and 0.0109 ± 0.0001 mm −1 in the in vivo forearm baseline measurements. Since the difference was negligible, the calibration equation derived using intralipid phantoms (figure 2) was employed for in vivo tissue temperature measurements. Finally, the corrected water spectra were fit to a library of pure water spectra at different temperatures to determine the temperature of the measured volume ( figure 3(b)). At this step, we found that the library should have pure water spectra at least at every degree Celsius for accurate temperature determination. The discrepancy observed at 750 nm between the tissue water spectrum and pure water spectrum in figure 3 might be due to over-correction of lipid and deoxy-hemoglobin in the tissue water spectrum. However, this wavelength range does not affect the temperature measurement since only wavelengths above 935 nm were employed for temperature calculation.
The algorithm for calculating absolute temperature was tested and validated in several different intralipid phantoms during heating and cooling using simultaneous DOS and thermistor measurements (figure 4). Excellent correlation was observed (R = 0.96) with a slope near unity (0.94) indicating that the sensitivity of optical recordings is comparable to the thermistor. The average difference in the measured absolute temperature between the two techniques was 1.1 ± 0.91 • C. This difference was similar to the standard error of the broadband DOS system (1 ± 0.74 • C). Thus, the difference between the two systems might be due to instabilities of the DOS system, likely because the upper usable limit of our spectrometer was about 1000 nm. As a result, in some cases the spectrum was shortened to less than 998 nm after the blue-shift and we extended it using an extrapolation function in the MATLAB. Because the sensitivity of the spectrometer in this wavelength region is poor, the extrapolated portion can have greater error. Spectra with very large errors in the extrapolation were excluded using an exclusion criteria based on approximate second derivatives ('diff' function of the MATLAB) of the range 981-998 nm. Any spectra with second derivatives greater than 7 × 10 −5 were excluded; however, many spectra with smaller but noticeable errors were included. When we used a different spectrometer with an upper limit of ∼1010 nm for the in vivo measurement, this problem did not appear. Thus, a spectrometer with the highest possible spectral range is recommended for temperature determination.
Despite these factors, overall accuracy was ∼1 • C regardless of whether we were monitoring heating or cooling. The scattering coefficients of our liquid phantom data were different during heating and cooling (for example, power law fit to μ s versus wavelength: y = 16648x −1.225 for heating and y = 20485x −1.285 for cooling, data not shown), demonstrating that irreversible structural/macromolecular changes in the phantoms during heating did not impact accuracy. This is due to the fact that temperature was determined directly from scattercorrected absorption spectra and therefore unaffected by scattering property changes. Other optical methods that are not model-based and do not correct for scattering properties may be subject to hysteresis effects during heating/cooling cycles (Belashev andTernovoi 1996, Kakuta et al 2008).
In general, the accuracy and dynamic range of these DOS temperature measurements can be improved by including more precisely measured pure water extinction coefficient spectra in the library of water spectra at various temperatures and by expanding the pure water spectral library to include a broader temperature range and finer thermal resolution. The fact that this method relies on the quality and number of water spectra in the library is a potential weakness. However, since the peak wavelength of the water spectrum decreases linearly as the temperature increases (R = −0.95, 95% confidence interval: −0.98 and −0.91, p-value <2.2 × 10 −16 in our water spectra library), we expect that if more precisely measured water spectra of higher and lower temperature are included in the library, both accuracy and dynamic range can be improved.
During in vivo tissue measurements, the initial temperature was about 28 • C, which is much lower than core body temperature (figure 5). This discrepancy is likely due to contributions from superficial tissues (i.e. skin and subcutaneous lipid) with thicknesses ranging from 2 to 5 mm in our subject population (Ducharme et al 1991). Forearm skin temperature has been measured in the range of 28.3-32 • C depending on location (Montgomery andWilliams 1976, Ducharme et al 1991). Although the forearm was immersed in 37 • C water for 5 min, the effect on tissue temperature was likely minor since approximately 20 min is required to achieve thermal equilibrium (Montgomery andWilliams 1976, Ducharme andTikuisis 1992).
The lowest tissue temperature measured (∼20 • C) after the 30 min immersion in cold water agreed well with previously published studies that report ∼20 • C measurements from an implanted multicouple probe in a forearm (Ducharme and Tikuisis 1992). When the cold water was replaced with the 37 • C bath, arm temperature eventually returned to baseline levels. The time required for the temperature decrease was about the same as recovery to the baseline. Although the arm was immersed for about 30 min in 37 • C water, the arm deep tissue temperature did not reach 37 • C. Thus, while the current DOS probe geometry (20 mm source-detector separation) has been shown to have the ∼6-10 mm mean interrogation depth, the temperature sensitivity map is not precisely characterized and could be weighted to more superficial structures.
The baseline and recovered temperatures were similar in all four subjects before and after the mild cold stress (27.6 ± 0.5 • C and 27.2 ± 0.3 • C, respectively). In three of the subjects, the temperature drop weakly correlated with an oxy-hemoglobin reduction (R = 0.55, 0.57 and 0.39, respectively in A, C and D). Oxy-hemoglobin changes are closely related to variations in blood delivery. Thus, the observed deep tissue temperature reductions during the cold stress may be accompanied by a decrease in vascular perfusion. In all subjects, the magnitude of the deoxy-hemoglobin change was smaller than the oxy-hemoglobin reduction (25.4 ± 8.7% compared to 52.9 ± 11.7%) suggesting that the rate of tissue oxygen utilization also decreases during cold stress although not much as the blood supply. During recovery, the deoxy-and oxy-hemoglobin concentration increased 36.9 ± 7.9% and 50.3 ± 15.8%, respectively, indicating re-perfusion and slightly increased metabolism.
In conclusion, we have developed a non-invasive optical method for absolute temperature measurements in thick tissues based on resolving opposing effects of water vibrational frequency shifts due to macromolecular binding. Unlike most reflectance methods for assessing temperature, this model-based approach accounts for thermally induced changes in light scattering and is essentially independent of denaturation of the test material. In tissue simulating intralipid phantoms, the difference between optical and conventional thermistor measurements was ∼1 • C and no hysteresis was observed between heating and cooling. In in vivo tissue measurements, temperature and hemodynamic variations were measured simultaneously in response to a mild cold stress. Because temperature, hemodynamic and metabolic parameters are co-localized, this non-invasive technique could potentially be employed to optimize dosimetry during thermotherapies and enhance the information content of metabolic diagnostics such as cold stress. In addition, because DOS is compatible with diffuse optical tomography (DOT), this technique could lead to the development of new thermal contrast elements for thick tissue functional tomography.