2.1.

Theory

The light-induced rate of heat energy $Q$ ($\mathrm{W}/{\mathrm{m}}^3$) generated in a tissue or material with absorption coefficient ${\mu }_{\mathrm{a}}$ equals $Q={\mu }_{\mathrm{a}}\mathrm{\Phi }$, where $\mathrm{\Phi }$ is the fluence rate. Temperature change caused by such a heat source is described by the heat diffusion equation:

Calculation of the time derivative $dT/dt$, and especially, the second derivative ${\nabla }^{2}T$ is challenging, due to the noise in the MRT data. The derivatives $dT/dt$ and ${\nabla }^{2}T$ were obtained separately by calculating the analytical derivative of local polynomial fits, performed at each pixel and each time point. To reduce computation time, the derivatives were calculated only within a spatial mask, defined by applying a threshold to the magnitude images of the thermometry data. For the time derivative at time point $t_i$, a second order polynomial fit was performed to a window of $N$ time points around $t_i$. For simulations and in vivo experiments, $N=15$ was used, while $N=11$ was used for phantom experiments. The time curves were split in three periods: before, during, and after laser irradiation. At the transitions between the periods, only the available time points from within that period were used for calculating the fits, because the curves are nondifferentiable at the transition. The polynomial fit functions and their analytical derivatives were evaluated at each time point, to obtain smoothed temperature curves and temperature derivative curves. The root mean squared error (RMSE) of the smoothed versus the original temperature curves was calculated, and a threshold was applied to the RMSE to exclude pixels with a low fit quality from the mask. Next, 3-D polynomials of the form $T=a_0+a_{1}x+a_{2}y+a_{3}z+a_{4}xy+a_{5}xz+a_{6}yz+a_7{x}^2+a_8{y}^2+a_9{z}^2$ were fitted to the smoothed temperature data within a spherical neighborhood of radius $R$ around each pixel. For simulations and in vivo experiments, $R=3.75\mathrm{mm}$ was used, while $R=3.0\mathrm{mm}$ was used for phantom experiments. To improve the signal-to-noise ratio and fit quality, pixels outside the mask were excluded from the neighborhood. The Laplacian at each point was obtained by calculating ${\nabla }^{2}T=2(a_7+a_8+a_9)$.

2.2.

Simulations

The feasibility of our method was tested on a 3-D virtual temperature dataset, obtained by simulating the heat generation and transport by a predefined fluence rate distribution in a nonscattering phantom with known thermal properties and absorption coefficient ${\mu }_{\mathrm{a}}$. A point light source, centered in a $30\times 30\times 30{\mathrm{mm}}^3$ volume of interest, was defined on a tetrahedral mesh of 80,412 points. Between $t=120$ and $t=420\mathrm{s}$ of the simulation, the light source was switched on with a fluence rate distribution of $\mathrm{\Phi }(r)=A\exp (-{\mu }_{\mathrm{a}}r)/r$, with amplitude $A=12.5\mathrm{W}/\mathrm{m}$ and $r$ the radial distance from the point source. In an infinite homogeneous medium with absorption coefficient ${\mu }_{\mathrm{a}}=0.5{\mathrm{cm}}^{-1}$, such a light source produces heat with a power density of $Q(r,t)={\mu }_a\mathrm{\Phi }(r)$ for $120\le t\le 420$, and $Q(r)=0$ otherwise. The heat transfer equation ($\rho c\frac{dT}{dt}=k{\nabla }^{2}T+Q$) was solved using finite-element method functions in the partial differential equation (PDE) toolbox of MATLAB R2016a. Temperature maps were generated from $t=0$ to 1020 s in steps of 10 s. In five slices, centered at $z=-2.2$, $-1.1$, 0, 1.1, and 2.2 mm around the point source, the temperature maps were interpolated to a Cartesian grid of $64\times 64\text{pixels}$, resulting in 0.47-mm spatial resolution, similar to our in vivo measurements. A second temperature dataset was obtained by adding Gaussian noise with a standard deviation of 0.43°C, to assess the effect of noise on the fluence rate calculations. Subsequently, the fluence rate calculation was performed on the noiseless and noisy temperature datasets, and the result was compared to the input fluence rate.

2.3.

Fluence Rate Probe

A model IP85 isotropic light probe (Medlight, Ecublens, Switzerland) was used for validation of fluence rate estimations in phantom and in vivo experiments. The probe consists of a ceramic spherical bulb mounted onto the distal end of a bare optic fiber. A fraction of the light collected by the bulb was guided through the fiber and detected by a photodiode on the proximal fiber end. The bulb collected light from all incident angles with nearly identical efficiency, except for the blind spot created by the fiber. The overall efficiency of the probe was calibrated, by placing the fiber in a light beam with known irradiance. Since the efficiency of the probe depends on the refractive index of its ambient medium,¹³ the calibration was done in water, which has a similar refractive index as tissue. The light was produced by a fiber-connected 655-nm laser (model WSLB-650-002-H, WaveSpectrum, Beijing, China). The fiber output was placed in front of an optical diffuser, to generate a 20-deg diverging beam with uniform light intensity profile. The irradiance of the beam was first quantified using a photodiode (model S121C, Thorlabs, Newton, New Jersey). Next, the photodiode was removed and the beam was aimed orthogonally onto the surface of a degassed water bath. The fiber was placed horizontally just below the water surface, with the bulb in the center of the beam, at the same longitudinal distance as the photodiode. The water container was made of black plastic, to minimize light reflections. The assumed irradiance at the probe was corrected for 2% reflectivity at the air–water surface. Since the probe was only slightly below the water surface and the divergence of the beam was small, no correction was applied for the change in beam diameter due to refraction. The optical power at the proximal end of the probe fiber was measured with a photodiode (model S150C, Thorlabs). A calibration factor was obtained by dividing the irradiance at the probe by the measured fiber power output.

2.4.

Phantom Preparation

Gel phantoms were prepared by adding 4% w/v gelatin from porcine skin (G2500, Sigma Aldrich, St. Louis, Missouri) to ultrapure water at 30°C while continuously stirring. Desired amounts of Evans Blue (E2129, Sigma Aldrich) and anatase titanium dioxide (232033, Sigma Aldrich) were added, as absorbing and scattering agents, respectively. Next, 0.1% v/v formaldehyde (252549, Sigma Aldrich) was added to facilitate cross-linking and the solution was heated to 45°C. It was then poured into the desired container or mold, which was placed in a vacuum chamber at 150 mbar for 15 min to remove air bubbles. The samples were stored in a 4°C refrigerator to fully crosslink overnight and used the next day.

2.5.

Phantom Characterization

The absorption coefficient as a function of Evans Blue concentration in the gelatin samples was determined by transmission spectroscopy. Samples were prepared in duplo, as described above in cuvettes (7590-05, Brand, Wertheim, Germany) with a path length of $d=10\mathrm{mm}$. Evans Blue concentrations of 2.53, 5.05, 10.1, 20.3, and $40.5\mathrm{mg}/\mathrm{L}$ were prepared, and no ${\mathrm{TiO}}_2$ was added. Two reference samples were prepared in the same way, but without Evans Blue. In a dark room, the cuvettes were placed in a cuvette holder (Ocean Optics, Dunedin, Florida), in between opposing ends of two optical fibers. One fiber was connected to a white light source (DH2000, Ocean Optics), while the other one was connected to a spectrophotometer (AvaSpec-ULS2048, Avantes, Apeldoorn, The Netherlands). Transmittance $T$ at 655 nm was calculated by dividing the signal intensities of each sample by that of a reference sample, after subtracting a dark measurement signal from both sample signals. The absorption coefficient in ${\mathrm{cm}}^{-1}$ was then calculated as ${\mu }_{\mathrm{a}}=-\ln T/d$.

2.6.

Phantom Magnetic Resonance Experiments

MRI experiments were performed with a 7 T scanner (BioSpec 70/30 USR, Bruker) using a quadrature 72-mm-diameter transmit and receive birdcage coil. $T_2$-weighted anatomical reference images were acquired using a multislice spin echo scan, with five axial slices centered on the light spot. Scan parameters were as follows: repetition time $(\mathrm{TR})=1000\mathrm{ms}$, echo time $(\mathrm{TE})=30\mathrm{ms}$, $\text{matrix}=128\times 128$, $\text{slice thickness}=1.0\mathrm{mm}$, $\text{slice gap}=0.1\mathrm{mm}$, field of view $(\mathrm{FoV})=6\times 6{\mathrm{cm}}^2$. The same geometry was used for all scans. $T_1$ and $T_2$ values of the phantom were determined with RARE relaxometry acquisitions, with the following parameters: RARE $\text{factor}=2$, $\mathrm{TE}=11$, 33, 55, 77, 99 ms, $\mathrm{TR}=593$, 813, 1096, 1491, 2155, 5000 ms. $T_2^{*}$ values were determined with a multigradient echo sequence, with the following parameters: ${\mathrm{TE}}_{\min }=2.9\mathrm{ms}$, $\mathrm{\Delta }\mathrm{TE}=3.5\mathrm{ms}$, ${\mathrm{TE}}_{\max }=60.4\mathrm{ms}$, flip angle $(\mathrm{FA})=30\mathrm{deg}$, $\mathrm{TR}=1500\mathrm{ms}$. For MRT, a FLASH sequence was used with TE approximately equal to $T_2^{*}$ (typically between 7 and 10 ms) and FA equal to the Ernst angle (typically between 12 and 14 deg). Other scan parameters were $\mathrm{TR}=76\mathrm{ms}$ and scan $\text{time}=9.7\mathrm{s}$. After acquisition of five baseline scans, the thermometry acquisition was repeated for the duration of the heating experiments.

A custom-built setup was used to perform laser-induced heating of the phantom inside the MRI scanner. The light from a 655-nm laser (model WSLB-650-002-H, WaveSpectrum, Beijing, China) was guided into the magnet by an optical fiber, collimated to an approximately parallel beam, and reflected downward by a prism onto the phantom. The diameter of the projected spot could be adjusted by changing the distance between fiber and collimator lens, and was between 0.65 and 0.90 cm in all experiments. Before each experiment, power output of the collimated beam was measured with an optical power meter. The spot size of the light bundle was measured with a ruler, and the laser current was adjusted to achieve the desired irradiance, assuming a uniform beam profile. Irradiance settings of 400 and $576\mathrm{mW}/{\mathrm{cm}}^2$ were used. During the heating experiments, the laser controller (model ITC4005, Thorlabs) was switched on after 1 min and off after 11 min, and was triggered by the MR pulse sequence. The complete experiment duration was 31 min, acquiring thermometry maps during a 20-min cooling period as well.

$B_0$ magnetic field drift correction was performed based on a linear space dependent phase correction map from the MR signal phase images, similar to De Poorter et al.¹⁴ Four horizontally positioned cuvettes filled with mineral oil were positioned around the gel sample in the corners of the imaging FoV. The phase correction map was obtained by a least squares fit plane through the phase of the four reference phantoms. The phase in the phantom was corrected at each time point by subtracting the phase correction map. Temperature maps were obtained by calculating the phase difference images $\mathrm{\Delta }\phi$ with respect to the average of the five baseline scans (${\phi }_0$), and division by the PRF change coefficient $\alpha =0.01\mathrm{ppm}/°\mathrm{C}$, the gyromagnetic ratio $\gamma$, the echo time TE, and the magnetic field strength $B_0$, according to $\mathrm{\Delta }T=\mathrm{\Delta }\phi /(\alpha \gamma B_0\mathrm{TE})$.

2.7.

Phantom Validation

For validation of MRT-based fluence rate results, fluence rate probe measurements were performed in similar phantoms as used in the MR experiments. The same optical setup was used to produce similar beam geometry. For quantification of the fluence rate along lines perpendicular to the optical axis of the light beam, the probe was inserted 2 cm into the side of the phantom, while the beam was projected onto the top surface. Using a motorized linear stage, the optical setup was then moved over the phantom in steps of 0.5 mm, such that the beam center traveled a straight path right above the probe tip. The depth $z$ of the probe with respect to the phantom’s top surface was measured with a caliper after slicing the phantom near the probe tip. The experiment was repeated in multiple phantoms, each time with the probe at a different depth. For measurement of the fluence rate profile along the optical axis, the probe was inserted into the bottom of the phantom, while the beam was again aimed onto the top surface and centered above the probe. The probe was inserted further into the phantom, in steps of 0.5 mm, until the tip emerged from the top surface.

2.8.

Mouse Cancer Model

All animal experiments were approved by the Animal Care and Use Committee of Maastricht University (protocols: 2012-139 and 2014-011). The maintenance and care of the experimental animals were in compliance with the guidelines set by the institutional animal care committee, accredited by the National Department of Health.

CT26.WT murine colon carcinoma cells (American Type Culture Collection (ATCC), CRL-2638) were cultured at 5% ${\mathrm{CO}}_2$ and 37°C in RPMI-1640 medium (Invitrogen, Breda, The Netherlands), supplemented with $50\mathrm{U}/\mathrm{ml}$ penicillin/streptomycin (Lonza Bioscience, Basel, Switzerland) and 10% fetal bovine serum (Greiner Bio-One, Alphen a/d Rijn, The Netherlands). Early passages (9-13) of the original ATCC batch were used for inoculation and $2\times {10}^6$ cells were inoculated subcutaneously in the right hind limb of 10- to 12-week-old BALB/c mice (Charles River, Maastricht, The Netherlands). Tumors became palpable at 3 to 5 days after inoculation. Experiments were performed when tumor volumes were between 271 and $621{\mathrm{mm}}^3$, approximately after 10 days. During MR imaging, mice were anesthetized using isoflurane (3.5% for induction, 1% to 2% for maintenance) in medical air flowing at $0.4\mathrm{L}/\min$. Breathing rate was monitored with a pressure balloon. Body temperature was monitored using a rectal probe, and kept at 36°C to 37°C using a warm water circuit. In the course of the study, $n=9$ mice were used, of which some were scanned two times on different days.

2.9.

In Vivo Experiments

The tumor-bearing hind limb was shaved before the experiment. A 19G needle was used to make a small perforation in the skin to facilitate insertion of a probe. Mice were placed in the same setup as used for phantom experiments. The hind limb was gently fixed in a stretched position, with the tumor on top and as close as possible to the isocenter of the MR scanner. An isotropic fluence rate probe was inserted in the tumor through the skin perforation. The penetration depth of the probe was deduced from the position of a pen marking on the fiber. Similar to the phantom experiments, a beam of 655-nm light was projected downward onto the tumor. The center of the laser beam was positioned in the same axial slice as the probe. Irradiances of 300, 400, and $500\mathrm{mW}/{\mathrm{cm}}^2$ were applied for 5 min, and the beam diameter was between 0.65 and 0.90 cm. The central MR slice was longitudinally aligned with the tip of the probe. An absorption coefficient of ${\mu }_{\mathrm{a}}=1.0{\mathrm{cm}}^{-1}$ was assumed for the fluence rate calculation, which is within the range of values for tumors reported in literature.^12,15

In order to demonstrate the MRT technique, a custom fiber optic temperature probe (Neoptix, Quebec, Canada) embedded in a 1-mm-diameter carbon needle was inserted into the tumor of one mouse, instead of a fluence rate probe. For this animal, the duration of laser heating was 13 min at an irradiance of $0.71\mathrm{W}/{\mathrm{cm}}^2$ and a beam diameter of 0.8 cm.

3.1.

Magnetic Resonance Thermometry Demonstration

Within the MRI research field, MRT is an established technique, although there are only a few demonstrations of MRT in small animals at 7 T.¹⁶ To assess accuracy and precision of the MRT, the MRI-based temperature data were compared to measurements with the tumor interstitial temperature probe (Fig. 1). The standard deviation of the temperature in a tumor region of interest (RoI) in a preheating temperature map was 0.43°C. Since MRT is a subtraction-based method, only temperature changes can be measured, relative to the starting temperature. Therefore, the robustness of the measurements decreases with the duration of the experiment, despite the use of field drift correction. For example, in Fig. 1(c), the RMSE of MRT versus the temperature probe increases from 0.49, to 0.90, and 1.29°C in the 1st, 15th, and 30th min, respectively. The temperature stability of in vivo experiments without heating was assessed by calculating the average absolute temperature drift in the final (20th) minute in a tumor RoI. The mean and standard deviation of this temperature stability over all in vivo experiments included in this study ($n=10$) was $1.20\pm 0.76°\mathrm{C}$.

Fig. 1

Demonstration of the MRT technique. (a) Anatomical reference image of the axial cross-section of a mouse hind limb. The red contour indicates the outline of the tumor. (b) Example of a temperature map at $t=14\min$, showing the temperature increase after 13 min of laser-induced heating. The temperature map shows a strong increase in the top of the tumor, where the light entered the tissue. The white arrow indicates the position of the temperature probe. (c) In a pixel next to the probe, the MRT data followed the probe temperature with an RMSE of 0.74°C, calculated over the entire duration of the experiment. The temperature starting point of the MRT data was set equal to the first measurement value of the probe.

3.2.

Simulations

To demonstrate the feasibility of MRT-based fluence rate calculation, 3-D simulations were first performed using virtual thermometry datasets, calculated for a known predefined fluence rate. Figure 2(a) provides a schematic outline of the simulation, with snapshot images of the central slices of the fluence rate distribution, temperature maps, the virtual MRT maps, and the final fluence rate map. The maximum temperature elevation at the center was 5.7°C. For the noiseless case, the spatial RMSE of the resulting fluence rate was only $8.8\mathrm{mW}/{\mathrm{cm}}^2$ for pixels located farther than 1.5 mm from the point source and calculated over the entire FoV, whereas the predefined fluence rate was larger than $100\mathrm{mW}/{\mathrm{cm}}^2$ for all pixels at a distance within 8.3 mm from the point source. The largest spatial RMSE was found close to the light source center, where the fluence rate asymptotically approaches infinity. Even with added Gaussian noise with a standard deviation similar to a typical in vivo MRT measurement (${\sigma }_{\mathrm{T}}=0.43°\mathrm{C}$), the fluence rate could be estimated with acceptable confidence (spatial $\mathrm{RMSE}=66.2\mathrm{mW}/{\mathrm{cm}}^2$). In Fig. 2(b), the time curve of the fluence rate in a pixel close to the virtual light source is shown (top), as well as a spatial profile along the radial axis of the light source (bottom). For the time curves, the temporal RMSE of the estimated versus the actual fluence rate was 21.8 and $60.6\mathrm{mW}/{\mathrm{cm}}^2$, for the simulation without and with added noise, respectively. For the spatial profile, the spatial RMSE’s were 24.4 and $36.7\mathrm{mW}/{\mathrm{cm}}^2$, for the noiseless case and the added-noise case, respectively.

Fig. 2

(a) Flow chart of simulation and example data showing the feasibility of the fluence rate calculation. A finite element method was used to solve the heat diffusion partial differential equation, to generate a spatial–temporal dataset of the temperature increase induced by an analytically described light source. The temperature dataset was solved on a fine tetrahedral mesh and then resampled to a Cartesian grid. By adding Gaussian noise ($\sigma =0.43$) to this dataset, a virtual MRT dataset was obtained. Fluence rate calculation was performed on this dataset, resulting in an estimated light-source map. The example shows a good resemblance with the ground truth fluence rate, even with added Gaussian noise. (b) In the top figure, the ground truth and estimated fluence rate in a single pixel are shown as a function of time, based on MRT data with and without added noise. In the bottom, the fluence rate profiles along a line through the point source are shown.

3.3.

Phantom Experiments

To calculate the dependence of the absorption coefficient ${\mu }_{\mathrm{a}}$ on the Evans Blue concentration in 4% gelatin phantoms, a transmission spectroscopy experiment was performed. The absorption coefficient at 655-nm wavelength increased linearly with the concentration of Evans Blue (Pearson’s $\rho =0.9998$). A linear fit through all data points led to the following equation: ${\mu }_{\mathrm{a}}=94.2[\mathrm{EB}]+0.07{\mathrm{cm}}^{-1}$, with [EB] the concentration of Evans Blue in $\mathrm{g}/\mathrm{L}$.

Next, the fluence rate estimation method was tested in a nonscattering gelatin phantom, which was heated by a laser beam while temperature changes were measured by MRT. The result is displayed in Fig. 3, displaying the calculation steps from temperature map to a fluence rate map. The temperature strongly increased where the light entered the phantom, but over time the rate of increase declined [Figs. 3(a) and 3(e)]. When the laser was switched on, the time derivative was positive and gradually approaching zero [Figs. 3(b) and 3(f)]. After switching off the laser, the time derivative became negative as the sample cooled down again. The Laplacian [Figs. 3(c) and 3(g)] was negative in most of the irradiated region during illumination but quickly went back to zero after the laser was turned off. The calculated fluence rate map [Fig. 3(d)] displayed a positive band along the beam axis, with decreasing values as function of penetration depth. The calculated fluence rate in a pixel near the probe [Fig. 3(h)] was close to the probe value during laser irradiation and returned to practically zero afterward. The median fluence rate in this pixel was $266\mathrm{mW}/{\mathrm{cm}}^2$, which was an underestimation of 27% compared to the median probe value of $363\mathrm{mW}/{\mathrm{cm}}^2$.

Fig. 3

Demonstration of the fluence rate calculation in a gelatin phantom with Evans Blue and no scattering material. The top row contains maps of (a) the temperature increase $\mathrm{\Delta }T$, (b) its time derivative $dT/dt$, (c) the Laplacian ${\nabla }^{2}T$, and (d) the resulting fluence rate $\mathrm{\Phi }$, at a time point just after the laser was turned on. The bottom row (e–h) shows time curves of the same quantities, in the pixel indicated by the red crosshairs in the $\mathrm{\Delta }T$ map. The red line in the time plots indicates the time point to which the maps correspond. The color bar applies to all four maps, and the range of each map is the same as the vertical axis of the corresponding time curve. The gray box indicates the time interval during which the laser was turned on. The black curve in the fluence rate time plot shows the fluence rate measured by the probe. The gap in the maps reveals the location of the probe.

Further validation of the MRT-based fluence rate estimation was performed in a nonscattering gelatin phantom ($20.3\mathrm{mg}/\mathrm{L}$ Evans Blue, ${\mu }_{\mathrm{a}}=1.98{\mathrm{cm}}^{-1}$) and a turbid gelatin phantom ($5.04\mathrm{mg}/\mathrm{L}$ Evans Blue and $2.50\mathrm{g}/\mathrm{L}$ ${\mathrm{TiO}}_2$, ${\mu }_{\mathrm{a}}=0.544{\mathrm{cm}}^{-1}$). Both phantoms were irradiated for 5 min with a 0.8-cm wide collimated beam, normally incident onto the top surface. Figure 4(a) shows the fluence rate map for the nonscattering phantom. In Fig. 4(b), the fluence rate as a function of penetration depth $z$ along the beam axis was compared to the theoretically expected profile based on Beer–Lambert’s law: $\mathrm{\Phi }(z)={\mathrm{\Phi }}_0\exp (-{\mu }_{\mathrm{a}}z)$. ${\mathrm{\Phi }}_0$ is the fluence rate value in the top layer ($z=0$) of the phantom, and was obtained by reducing the irradiance ($E=400\mathrm{mW}/{\mathrm{cm}}^2$) by the specular reflection loss, expressed by the reflectance $R_{\mathrm{sp}}={[(n_0-n_1)/(n_0+n_1)]}^2$, where $n_0$ and $n_1$ are the refractive indices of air and the gelatin mixture, respectively. Assuming that $n_0=1.0$ and $n_1=1.33$, 98% of the irradiance was expected to be transmitted into the phantom, so ${\mathrm{\Phi }}_0=392\mathrm{mW}/{\mathrm{cm}}^2$. Based on the Evans Blue concentration of the phantom, the expected absorption coefficient was ${\mu }_{\mathrm{a}}=1.98{\mathrm{cm}}^{-1}$. The RMSE of the experimental versus the theoretical fluence rate along the beam axis was $23.6\mathrm{mW}/{\mathrm{cm}}^2$ (6.0% compared to ${\mathrm{\Phi }}_0$), and $R^2$ was 0.929. The experimental fluence rate generally was a bit lower than theoretically predicted, especially in the top layer of the phantom.

Fig. 4

Validation of MRT-based fluence rate maps of two laser-irradiated phantoms. (a) Fluence rate map for a gelatin phantom with $20.3\mathrm{mg}/\mathrm{L}$ Evans Blue and without a scattering material. (b) Comparison of the fluence rate attenuation along the optical axis (blue solid line) to the expected exponential decay profile (red dashed line). (c) Fluence rate map for a gelatin phantom containing $5.04\mathrm{mg}/\mathrm{L}$ Evans Blue and $2.5\mathrm{g}/\mathrm{L}$ ${\mathrm{TiO}}_2$. (d) Comparison of the fluence rate attenuation along the optical axis (blue solid line) to probe measurements (red dashed line). (e) Comparison of the MRT-based fluence rate (solid lines) to probe measurements (dashed lines) along radial lines at four different depths $z$ below the phantom surface. The two fluence rate maps show the median of the fluence rate over the interval during which the laser was turned on. Note that the color bars of both phantoms have different ranges.

The fluence rate map of the scattering phantom [Fig. 4(c)] shows a shallower light penetration than for the nonscattering phantom. The fluence rate in the top layer exceeded the applied irradiance values. Figures 4(d) and 4(e) show comparisons between the MRT-based fluence rate and probe measurements along the beam axis and along radial lines at four different depths. The spatial RMSE for the longitudinal profile was $51.1\mathrm{mW}/{\mathrm{cm}}^2$ (7.0%), whereas for the radial profiles, it was 129.4 (17.7%), 34.5 (4.7%), 35.7 (4.9%), and 34.7 $(4.7\%)\mathrm{mW}/{\mathrm{cm}}^2$ at $z=1.0$, 2.0, 4.8, and 7.2 mm, respectively. The percentages are all with respect to the maximum probe fluence rate along the optical axis. $R^2$ values were 0.942 for the longitudinal profile, and 0.714 and 0.918 for the radial profiles $z=1.0$ and 2.0 mm. For the other two depths, $R^2$ was not calculated, because the fluence rates were approximately zero.

3.4.

In Vivo Experiments

Finally, the method was demonstrated in vivo in a mouse tumor model. An example of these experiments is shown in Fig. 5. Essentially, the observations were similar to those in the phantom experiment. The overall fluence rate distribution in the tumor was as expected: high values in the upper region, where the beam entered the tissue, but a more diffuse profile than for the nonscattering phantom. In most pixels, the fluence rate time curves were close to zero when the laser was off, and increased when it was on. In general, the temperature SNR of the in vivo experiments was lower than the phantom experiments. Signal voids were observed in some tumors, especially when they were scanned for the second time, possibly due to hemorrhages induced by the probe or laser heating. In addition, some temperature drifts were observed after laser irradiation in regions of most tumors, e.g., slow increases after initial decrease or decrease below zero. One mouse was excluded from analysis, because the signal void around the probe was too large to reliably calculate ${\nabla }^{2}T$ in large parts of the tumor. The correlation between the fluence rate measured by the probe and the MRT-based fluence rate near the probe was statistically significant ($\rho =0.896$, $p=0.0026$), as shown in Fig. 6. The median value of the fluence rate over the interval during laser irradiation was used. A linear fit through the data points resulted in a slope somewhat below unity, which was caused mostly by the two data points with the largest probe values. All the other data points were closer to the identity line.

Fig. 5

Example of in vivo fluence rate calculation. The top row shows maps of (a) the temperature increase $\mathrm{\Delta }T$, (b) its time derivative $dT/dt$, (c) the Laplacian ${\nabla }^{2}T$, and (d) the resulting fluence rate $\mathrm{\Phi }$, at a time point just after the laser was turned on. The maps are presented as color overlays over a grayscale anatomical reference image. The bottom row (e–h) shows time curves of the same quantities, in the pixel indicated by the red crosshairs in the $\mathrm{\Delta }T$ map. The probe was located just below this pixel. The red line in the time plots indicates the time point to which the maps correspond. The colorbar applies to all four maps, and the range of each map is the same as the vertical axis of the corresponding time curve. The gray box indicates the time interval during which the laser was turned on. The black curve in the fluence rate time plot shows the fluence rate measured by the probe.

Fig. 6

Scatter plot of the MRT-derived fluence rates of in vivo experiments, compared to the probe values. Each dot represents one experiment. The correlation coefficient is displayed above the figure and was statistically significant. The blue solid line is a linear fit through the origin. The identity line is displayed as a dashed black line.