1.

Introduction

Diffuse optical tomography (DOT) is an emerging functional imaging technique with a great potential in medical imaging.^1–5 It uses near-infrared (NIR) light to probe deep tissues because of the low-optical attenuation in this spectral window.⁶ By imaging metabolic compounds such as water, fat, and oxy- and deoxy-hemoglobin,⁷ DOT is an emerging imaging tool for breast cancer monitoring^8–19 and functional brain imaging.^20–25 Despite extensive efforts, its translation to the clinical arena has been delayed due to its poor spatial resolution. In fact, the highly scattering nature of the biological tissue makes the reconstruction of DOT images a very challenging process. Several factors degrade the spatial resolution and quantitative accuracy of DOT. Some of these are the ill-posed nature of the inverse problem and nonuniqueness of the solution, which results from the fact that data are measured only at the boundary of the imaged medium.²⁶ Also, solving the inverse problem in DOT requires the inversion of the sensitivity matrix describing the relationship between each measurement and the unknowns to be reconstructed. Since the number of measurements is much smaller than the unknowns, DOT inverse problem is underdetermined and the sensitivity matrix is nonsquare.²⁷ Therefore, its inversion is not a straight forward process and requires the use of unstable pseudo inversion methods.^28,29 It has been demonstrated that the quality of reconstructed tomographic images can be considerably improved in general, if internal measurements are available.³⁰ In fact, this reduces the under-determination of the inverse problem and makes it well-posed. Accordingly, the ideal scheme for DOT would be performing optical measurements from the whole medium under investigation, which is currently not possible with conventional DOT approaches.

Meanwhile, an intriguing combination of optical and ultrasound techniques has led the development of photoacoustic imaging (PAI) that can provide optical absorption maps with much higher resolution with several centimeters of depth penetration.^31–33 In PAI, a short-pulsed laser ($\sim 10\mathrm{ns}$) is used to irradiate the sample. The sound wave produced by the thermoelastic expansion of the tissue as a result of the energy deposition of the incident laser light is detected by ultrasound transducers. The ultrasound waves scatter much less in biological tissue than optical waves and this drastically improves the overall resolution. In addition to requirement of a pulsed laser, PAI requires transducer-tissue contact due to utilization of sound waves. Although great deal of progress has been achieved in this powerful imaging technology, thick tissue imaging is challenging and considerable effort is still being spent for its clinical translation.³⁴

Previously, we introduced a new alternative technique that we termed photomagnetic imaging (PMI).^35,36 Although PMI aims to recover spatially resolved optical absorption of the medium under investigation, it does not utilize any optical detectors. Instead of the conventional boundary optical measurements, PMI utilizes magnetic resonance thermometry (MRT) to measure internal spatiotemporal distribution of temperature variation induced by the local absorption of light when the medium is illuminated with a laser. Despite its similarity to PAI in terms of leveraging temperature modulation induced by laser illumination, PMI is inherently different due to utilization of noncontact MRT measurements and slow heating process achieved in a much longer time-frame (i.e., tens of seconds). PMI provides high-resolution optical absorption maps using the measured temperature maps and modeling of light propagation in addition to heat diffusion in tissue. The dedicated PMI reconstruction algorithm is based on the minimization of the difference between the measured and simulated spatiotemporal temperature maps. These simulated time-dependent temperature maps are obtained by solving the coupled diffusion and Pennes bioheat equations with finite element method (FEM).^2,37–39 In fact, this system of equations models the photon density distribution in the medium and the conversion of photons energy to bioheat by considering the local optical absorption.^36,40

Since the quality of the reconstructed PMI optical absorption images directly depends on the quality of the temperature measurements, the choice of MRT imaging parameters is critical. Several MRT parameters can be used to determine the absolute or relative temperature inside a medium such as T1 or T2 relaxation time, the diffusion coefficient, magnetization transfer, or proton resonance frequency (PRF). Unlike T1 or T2 MRT, PRF is linearly dependent on temperature. Moreover, PRF has been widely used to monitor the temperature variation when performing ablation therapy.⁴¹ Technically, the temperature increase causes the decrease of hydrogen bonding which decreases the PRF. Using any phase sensitive magnetic resonance imaging (MRI) sequence such as gradient echo, asymmetric spin echo, or Echo-planar imaging, a proportional relationship is established between the phase accumulation and temperature change between successive frames.^42,43 PRF is not only fast, realizing up to $10\text{frames}/\mathrm{s}$ but also very sensitive and invariant for different types of tissues.⁴⁴ Also, by reducing the temporal resolution, its sensitivity can be enhanced further.

In our previous work, we presented the feasibility of PMI.³⁵ In these initial studies, the measured temperature maps were matched with the estimated ones using the FEM-based multiphysics forward solver. This first prototype PMI system had a single illumination port and a low signal-to-noise ratio (SNR) that resulted in utilization of high-illumination power, nearly three times of the maximum skin exposure limits defined by the American National Standards Institute (ANSI). In this study, we present a full four-port illumination PMI system with a higher SNR which brings down the laser power levels below the ANSI limits and paves the way for the future clinical use of PMI for diagnostic purposes. More importantly, we present the first experimental studies that involve the reconstruction of optical absorption maps using a dedicated PMI reconstruction algorithm, which was tested earlier with simulation studies.³⁶ To evaluate the performance of the PMI in resolving inclusions as a function of their absorption contrast, size, and depth, multiple phantoms are constructed and utilized. The results confirmed that PMI is not only able to resolve inclusions as small as 1-mm diameter but also recover their absorption coefficients accurately.

2.

Method

3.

Experimental Studies

The MRT measurements are performed inside a Philips 3 Tesla Achieva system. PRF images are acquired using a gradient echo sequence using 60 and 12 ms as repetition and ET, respectively. Figure 1(a) shows a schematic of our special designed animal MRI coil inside the MR bore. The illumination of the phantom is realized using four 7-W laser diodes emitting at 808 nm (Focuslight, China). The laser diodes and their drivers are all positioned at the control room, far away from the magnet. Light is transported to the PMI interface located inside the MR bore using four 15-m long optical fibers. The PMI interface consists of a specially designed RF coil with four windows for illumination and four ports that hold the collimation optics (Fig. 1). When optical fibers are connected to these ports, their light output is collimated using 35-mm Newport optics aspherical lenses prior to reaching the sample positioned at the center of the RF coil. The laser power per unit area is set to the ANSI limits ($0.32\mathrm{W}/{\mathrm{cm}}^2$ for 808 nm) and confirmed by a power meter (PM100D, Thorlabs).

Fig. 1

(a) A schematic of PMI setup showing the phantom and the optical instrumentation inside the MRI bore. (b) The picture of the PMI interface sitting on the MRI bed. It consists of a specially designed RF coil with four windows for illumination and four ports that hold the collimation optics.

The MR temperature maps acquisition is synchronized with the laser drivers. Figure 2(a) shows the timeline of data acquisition. First, a T1 weighted low-resolution MR pilot image is acquired to localize the axial position of the laser probe. Once the axial plane located, a dynamic imaging set consisting of multiple frames (6 s each) is initiated. Gradient echo sequence is used to obtain a high-resolution phase image at each time point. The first baseline frame (${\phi }_0$) is acquired before turning-on the laser. After the lasers are turned on, the phantom is illuminated simultaneously from four sides during the heating phase, and another frame is acquired (${\phi }_i$) (Fig. 2). Therefore the total PMI data acquisition time is nearly 12 s. After subtracting the baseline image (${\phi }_0$) from the image acquired during heating phase (${\phi }_i$), relative phase change for each pixel can be obtained and hence, a high-resolution image describing the temperature increase can be attained using Eq. (5). In fact, each pixel of this image serves as a measurement point to monitor the temperature increase. This temperature map is utilized as the final measurement map in the PMI reconstruction minimization process [Eq. (1)]. Actually, more than one frame can also be acquired during the heating phase to boost the temperature change and hence, SNR. Although, the best PMI reconstruction should be obtained using phase image acquired at the end of the heating cycle in general, combination of several frames can also be utilized.³⁶ Note that PMI data can be acquired continuously as long as the laser illumination is kept under ANSI limits.

Fig. 2

(a) Timeline of PMI data acquisition showing the laser status, sample temperature at any particular point and the MRT acquisition for different cycles. (b) An illustration for the mesh of the cross section of the phantom and the collimated laser illumination beams directed onto the sample from four sides.

4.

Results

4.1.

Forward Problem Validation

The first step in the experimental validation of our technique is the validation of the forward problem. For this purpose, a PMI temperature map is first acquired using a mice-sized 25-mm diameter cylindrical agarose homogeneous phantom with an optical absorption coefficient of $0.01{\mathrm{mm}}^{-1}$, to mimic biological tissue. A 5-mm inclusion is embedded inside the phantom, 6-mm deep under the illumination site in order to mimic heterogeneous medium having a higher absorption area [Fig. 3(a)]. The inclusion is placed a bit off-center and its optical absorption is set to be eight times higher than the background.

Fig. 3

(a) The phantom cross section showing the inclusion size and position. The laser used to heat the phantom from its top side is represented by the red arrow S. The temperature maps (b) measured using MRT and (c) simulated using the forward solver. (d) The profiles taken along the $y$-axis on the measured and simulated temperature maps.

Meanwhile, using the forward solver, the temperature map is simulated on an identical synthetic phantom having the same size, optical, and thermal properties. This simulated temperature map is then compared to the temperature map obtained by MRT.

Figures 3(b) and 3(c) show the measured and simulated temperature maps, respectively. Note that, in order to compare the measured and the simulated data, a calibration step is needed. During this step, two parameters are defined from the experimental setup: the power of the laser and the shape of the illumination beam. First, the laser power is measured at the output of each of the four collimation lenses using a power meter. Second, PMI measurements are performed on a homogeneous phantom with known optical and thermal properties. The PMI interface holds the phantom at the center and therefore, the position of the illumination beams on the sample is approximately known. Nevertheless, a minimization program calculates the exact position of the illumination beams on the sample by minimizing a quadratic difference between the measured and simulated temperature maps. Using these calculated parameters, a very good agreement is achieved between simulated and measured maps on this heterogeneous phantom as shown by the profiles shown in Fig. 3(d).

These first simulations show that the variation in the temperature decreases quickly with depth and goes below the sensitivity of our instrumentation. Additional simulations are performed to analyze the effect of simultaneous multiple point illumination on the depth penetration and stability of the inverse problem rather than sequential illumination from each point. As seen in Fig. 4, as the number of illumination beams increases, higher SNR is obtained at the center of the phantom. The temperature maps reveal that, four-port illumination gives the highest overall temperature increase in the whole phantom. The profiles in both $x$- and $y$-directions, clearly show that the temperature increase at the center of the phantom is the highest and more importantly above the noise level, only when four-port illumination is used [Fig. 4(b)]. Based on these results, we upgraded our PMI interface by adding four-port illumination and performed all the experiments with this configuration.

Fig. 4

Temperature maps simulated when: (a) one, (b) two, and (c) four lasers are used. The temperature profiles along the (d) $y$-axis and (e) $x$-axis. The noise level of the MRT sequence corresponds to 0.1°C (green highlighted). With only one laser illumination from the top of the phantom, the signals measured at the bottom half of the phantom is less than the noise level (red). When two-port illuminations are utilized, the measurements at the center of the phantom suffer from low SNR in the direction of illumination and mostly under the noise level in the orthogonal direction (blue). By increasing the illumination ports to four, this problem is solved and the SNR is increased in the overall phantom (green).

4.2.

Phantom Experiments

Three different phantoms are imaged to validate the ability of our technique in resolving inclusions with different depths, sizes, and contrasts. All experiments are realized using a mice-sized 25-mm diameter cylindrical agarose phantom mimicking a small animal. The optical absorption coefficient of this cylindrical phantom is set to $0.01{\mathrm{mm}}^{-1}$. In the following experiments, different inclusions are implanted in this phantom in order to model tumors located at different depths, with various sizes and contrasts. The optical absorption coefficients of these inclusions are set according to the requirements of each experiment and are presented bellow. All the following reconstructions are made using a fine mesh consisting of 17,586 triangular elements connected at 8926 nodes. This mesh is generated using a maximum element size set to 0.25 mm which corresponds to a quarter of the size of the smallest inclusion to be reconstructed.

4.2.1.

Depth phantom

In this experiment, we embedded three 2.5-mm diameter inclusions with the same absorption coefficient set to be four times higher than the absorption of the background. The inclusions are placed at 3, 7.5, and 12 mm under the surface of the phantom as shown in Fig. 5(a). Figure 5(b) shows the measured temperature map obtained using MRT. This map clearly shows the increase in temperature induced by the laser under the illumination sites. Note that the temperature does not increase similarly under the four lasers illumination sites due to different power output levels of each port as well as variation in spot size. In practical applications, the sample will not be perfectly circular and hence the spot size as well as the power density might be different at each illumination site. Due to this reason, we intentionally set the power and spot size for each port to be a bit different in these experiments in order to model a realistic case. As mentioned above in Sec. 4.1, PMI calibration measurement is first performed on a homogeneous phantom with known optical and thermal properties. This calibration measurement is utilized to find the differences at each illumination site and hence to figure out any mismatch between the experiment and the simulation parameters by minimizing the objective function [Eq. (1)].

Fig. 5

(a) Axial MRI of the cylindrical phantom. The size and position of the inclusions are shown using black circles since no absorption contrast can be seen on the MRI image. All inclusions have the same contrast set to four times more absorbent than the background. (b) The temperature map at the end of the heating cycle. (c) The PMI absorption reconstructed map. (d) The absorption (red) and the temperature (blue) profiles across the three inclusions along the red arrow shown on (a).

Even without any calibration, Fig. 5(b) shows that the temperature map already reveals the position of the three inclusions prior to any reconstruction process. In fact, due to their higher optical absorption, the temperature increase is higher inside the inclusions. However, strong depth dependence is observed even though the three inclusions have identical size and contrast. Figure 5(d) shows clearly the observed depth dependence on the profiles (blue) calculated on the temperature map, along the red arrow shown in Fig. 5(a).

In order to overcome the depth dependence, and obtain quantitative absorption maps, the PMI inverse problem is solved. Figure 5(c) shows the reconstructed absorption map when the measured MRT temperature map is used as the input for the solver. The three inclusions are resolved accurately and their absorption coefficients are recovered without any depth dependence. Indeed, the three reconstructed inclusions have an average error of 4% in diameter and only 7% in the recovered absorption coefficient. The obtained results are summarized in Table 1.

Table 1

Summary of PMI experiments performed with the depth phantom. All three objects located at different depths are recovered successfully with 4% error in size and 7% error in absorption coefficient.

	Inclusion 1		Inclusion 2		Inclusion 3
	Real	Recon	Real	Recon	Real	Recon
Diameter (mm)	2.5	2.6	2.5	2.4	2.5	2.6
${\mu }_a$ (${\mathrm{mm}}^{-1}$)	0.04	$0.039\pm 0.002$	0.04	$0.037\pm 0.002$	0.04	$0.038\pm 0.002$

4.2.2.

Size phantom

In this case, four inclusions with different sizes are embedded $\sim 6\text{-}\mathrm{mm}$ deep into the 25-mm diameter cylindrical agarose phantom. The inclusions are 1.0-, 1.5-, 2.0-, and 3.0-mm diameter as shown in the cross section of the phantom shown in Fig. 6(a). In order to evaluate the ability of PMI in resolving inclusions with different sizes, we first chose the same optical absorption coefficient for the four inclusions and set it to be four times higher than the absorption of the background. Second, all four inclusions are placed at the same distance from the surface of the phantom. This ensures that each inclusion is at the same distance from the edge under the illumination in order to limit the depth dependence of the heating results. As mentioned above, the variation in the illumination beam at the four sides is corrected using the calibration procedure. Figure 6(b) shows the high-resolution absorption map obtained after solving the PMI inverse problem. The inclusions are recovered successfully on the reconstructed absorption map and their location shows a very good agreement with the cross section of the phantom shown in Fig. 6(a).

Fig. 6

(a) Axial MRI of the cylindrical phantom. The size and position of the inclusions are shown in black circles since no absorption contrast can be seen on the MRI image. All inclusions have the same contrast set to be four times more absorbent than the background. (b) The PMI absorption reconstructed map. (c) The absorption profiles across the four inclusions along the red lines shown on (a).

Figure 6(c) shows the four reconstructed absorption profiles taken along the red lines across each inclusion as shown in Fig. 6(a). The profiles show clearly the accuracy in recovering the sizes and contrasts of the inclusions. Table 2 summarizes the results of the experiments using the size phantom.

Table 2

Summary of PMI experiments performed with the size phantom. All four inclusions having different sizes are recovered successfully. Note that the error in the recovered absorption of the smallest inclusion (1 mm) is high since it has approximately the same size as the MRT pixel size (0.8 mm).

	Inclusion 1		Inclusion 2		Inclusion 3		Inclusion 4
	Real	Recon	Real	Recon	Real	Recon	Real	Recon
Diameter (mm)	1	1.1	1.5	1.6	2	2.0	3	3.0
${\mu }_a$ (${\mathrm{mm}}^{-1}$)	0.04	$0.031\pm 0.001$	0.04	$0.0383\pm 0.001$	0.04	$0.0391\pm 0.001$	0.04	$0.039\pm 0.003$

As seen from Table 2, the size of the largest two objects, inclusions 2.0 and 3.0 mm in diameter, are recovered perfectly, while their absorption coefficients are recovered with only 2% error. The error in the recovered size increased to 7% for the 1.5-mm diameter inclusion, while the error in the recovered absorption is only 5%. Although the size of the smallest inclusion, 1 mm in diameter, is recovered with only 10% error, the quantitative accuracy in the recovered absorption coefficient is the lowest where the error reached 32%. In fact, since the pixel size of the MRT image is set to 0.8 mm, the higher error in the recovered absorption coefficient of this 1 mm small inclusion is well expected. However, it is important to notice that although the size of this inclusion is very close to the MRT pixel size, it is accurately localized. This problem can be easily overcome by decreasing the size of the MRT pixel size but in expense of PMI temporal resolution and computation time.

4.2.3.

Contrast resolution

In this last experiment, four 5-mm diameter inclusions are embedded in the 25-mm cylindrical agarose phantom at the same distance, $\sim 5\mathrm{mm}$, under its surface. Taking the same inclusion size and distance under the surface for all inclusions reduces the size- and depth-dependence errors and allows us to focus on the ability of PMI in resolving inclusions with different contrasts. The optical absorption of the inclusions are set to be 2, 4, 6, and 8 times higher than the background as shown in Fig. 7(a).

Fig. 7

(a) Axial MRI of the cylindrical phantom. The size and position of the inclusions are delimited with the black circles since no absorption contrast can be seen on the MRI image. The inclusions are 2, 4, 6, and 8 times more absorbent than the background. (b) The PMI absorption reconstructed map. (c) A graphic showing the linear proportional relationship between the reconstructed and the real absorption coefficient.

Figure 7(b) shows the reconstructed absorption map which represents a high resolution PMI image of the phantom. This map shows a perfect agreement with the real absorption map presented in the cross section of the phantom [Fig. 7(a)]. The linear proportional relationship between the real and reconstructed contrasts of the four inclusions is shown in Fig. 7(c).

The results of this experiment are summarized in Table 3. The sizes of all four objects are recovered with $<3\%$ error. Similarly, the errors in the recovered absorption coefficients are all $<6\%$. Figure 7(c) shows the linear proportional relationship between the reconstructed and the real absorption coefficient demonstrating the high performance of PMI in recovering all four contrasts. Table 3 summarizes the results of the experiments using the contrast phantom.

Table 3

Summary of PMI experiments performed with the contrast phantom. All four objects with different contrasts are recovered successfully.

	Inclusion 1		Inclusion 2		Inclusion 3		Inclusion 4
	Real	Recon	Real	Recon	Real	Recon	Real	Recon
Diameter (mm)	5	5.1	5	4.9	5	4.8	5	4.8
${\mu }_a$ (${\mathrm{mm}}^{-1}$)	0.02	$0.019\pm 0.001$	0.04	$0.040\pm 0.002$	0.06	$0.059\pm 0.002$	0.08	$0.075\pm 0.004$

5.

Conclusion

In summary, we experimentally validated the ability of PMI in providing high-resolution absorption images from internal temperature measurements. PMI combines the high sensitivity of optical imaging with the high-spatial resolution of MRI. We have previously introduced the PMI imaging technique and demonstrated its quantitative imaging ability with simulation studies.³⁶ Here, we not only validate the superior performance of PMI with experimental studies but more importantly, present a four-port PMI system that can operate under ANSI limits ($0.32\mathrm{W}/{\mathrm{cm}}^2$ for 808 nm). That is a critical milestone since our ultimate aim is to utilize this imaging technology for diagnostic purposes such as breast or head and neck cancer imaging. Meanwhile, immediate application of the PMI technology can be preclinical imaging. For this reason, in this study, we choose to use mice-sized 25-mm diameter agar phantoms mimicking small animals.

Since we would like to keep the laser power level under ANSI limits, the challenging part of PMI is the sensitivity of the MRT. We found that the noise level of our MRT sequence is 0.1°C but we believe that it can be improved further by optimizing the MRT pulse sequence. Our experimental studies show that this noise level is enough to perform PMI imaging successfully.

Note that the SNR of the MRT signal from each pixel is depth independent. Meanwhile, the key measurement in PMI is the local temperature increase. Since the photon density decreases drastically with depth in turbid media, PMI probing depth directly depends on the sensitivity of the MRT. Therefore, the key strength of PMI is that as long as the temperature change is detectable by MRT, the resolution is preserved at any depth. Supporting this idea, our experimental results showed that PMI can recover 3-mm diameter inclusions buried at different depths successfully with only 7% error in the recovered absorption coefficient. The size phantom experiment also revealed the superior performance of PMI. Since MRT can detect the 1-mm inclusion located 6-mm deep, PMI can recover it successfully with only 32% error in the absorption coefficient.

Our experiments showed that we can probe 25-mm diameter sample successfully using four-port illumination. We believe that the probing depth can be increased using a higher number of illumination ports and ultimately developing a different scheme such as using a light guide surrounding the sample allowing an all-around homogeneous illumination. Meanwhile, we are working on improving the MRT parameters to increase its sensitivity and hence the probing depth of PMI. In addition, utilizing multiwavelength illumination, PMI can also reveal high-resolution images of tissue chromophores such as oxy- and deoxy-Hb, water, and fat. In fact, PMI can also be utilized with exogenous contrast agents such as gold nanoparticles and the Food and Drug Administration approved indocyanine green. Therefore, it will certainly make a high impact in the clinical management of cancer with the development of smart tumor targeting probes. For example, our preliminary simulation studies showed that this technique can be applied for breast cancer, if the breast is slightly compressed down to 4- or 5-cm thickness and illuminated from all around. We are currently performing extensive simulation and experimental studies to show the feasibility of this technique for breast cancer.