2.1.

Photon Migration in a Turbid Medium

Ultimately, the best way to provide a thorough description of a photon migration in a turbid medium is TD diffuse optical tomography, which decouples the underlying scattering and absorption of photons. NIR photon migration in a turbid medium is governed by the diffusion equation for the diffuse photon fluence rate $\phi (\mathbf{r},t)$ ,

2.2.

Light Propagation Model

A schematic diagram of the experimental setup of the eXplore Optix™-MX2 system (ART Advanced Research Technologies, Inc., Montreal, Canada) is shown in Fig. 1c . As shown in Fig. 1c, the whole light propagation procedure from source to detector can be divided into four steps: $\left[f_1\right]$ Green’s function representing a photon migration from source $\left({\mathbf{r}}_s\right)$ to a fluorescent inclusion $\left(\mathbf{r}\right)$ , $\left[f_2\right]$ lifetime decay $\left(\tau \right)$ of a fluorescent inclusion, $\left[f_3\right]$ another Green’s function showing a migration from a fluorescent inclusion $\left(\mathbf{r}\right)$ to detector $\left({\mathbf{r}}_d\right)$ , and $\left[f_{\mathrm{IRF}}\right]$ impulse response function (IRF) of approximated Gaussian wave form. In general, ${\mu }_a$ and ${\mu }_s^{\prime }$ are functions of wavelength, thus the steps $\left[f_1\right]$ and $\left[f_3\right]$ should be described by different propagators. However, since the excitation and emission wavelength are close to each other where tissue optical properties vary slowly, these Green’s functions can be described by the same propagator. A summary of the previous procedure can be written as the forms of the following analytical functions,

Fig. 1

(a) Reflection geometry is composed of a single source-detector pair ( $3\text{-}\mathrm{mm}$ separation), which is raster scanned every $2\mathrm{mm}$ across the phantom surface ( $3\times 3{\mathrm{cm}}^2$ area) encompassing the submerged fluorescent inclusion. (b) Side view of (a). (c) Schematic diagram of the experimental setup for the propagator description of photon migration.

To generate a fluorescent distribution in the region of interest, the integral in Eq. 4 is to be discretized [ $\mathbf{r}$ , $n\left(\mathbf{r}\right)$ , $\tau \left(\mathbf{r}\right)$ , ${\mu }_a\left(\mathbf{r}\right)$ , ${\mu }_s^{\prime }\left(\mathbf{r}\right)\to {\mathbf{r}}_j$ , $n_j\left(\mathbf{r}\right)$ , ${\tau }_j\left(\mathbf{r}\right)$ , ${\left[{\mu }_a\left(\mathbf{r}\right)\right]}_j$ , ${\left[{\mu }_s^{\prime }\left(\mathbf{r}\right)\right]}_j$ ]. For a series of measurements that have been made at source $\left[{\left({\mathbf{r}}_s\right)}_i\right]$ and detector $\left[{\left({\mathbf{r}}_d\right)}_i\right]$ positions, Eq. 4 turns to the following matrix equation,

3.1.

The eXplore Optix™-MX2 System

The eXplore Optix™-MX2 system is the only commercially available multiwavelength (470, 635, 670, and $760\mathrm{nm}$ ) TD optical molecular imaging system for small animals.¹³ The laser pulse width was approximately $70\mathrm{ps}$ and the temporal resolution of the detection system was approximately $250\mathrm{ps}$ . Schematic scanning geometry is depicted in Figs. 1a and 1b. Scanning resolution was $2\mathrm{mm}$ per step over the region of interest (ROI) with a scan time of $1\mathrm{s}$ per step. The ROI was a square of area $3\times 3{\mathrm{cm}}^2$ within which a small cylindrical Cy5 or Cy7 inclusion was submerged in the optical scattering Intralipid medium. The system uses a single source-detector configuration to scan the small animal in reflection mode, and measures the fluorescence TPSF for each scanned pixel. The fluorescence TPSFs were used to quantify $n$ , $\tau$ , and $d$ of the fluorescent inclusion.

3.2.

Phantom Manufacture

To expedite repeated and robust experiments, we manufactured solid optically scattering fluorescent pellets with either Cy5 or Cy7 at $10\mu \mathrm{M}$ , similar to a method that is found in the literature.¹⁴ The diameter and the height of these cylindrical pellets were 15 and $8\mathrm{mm}$ , respectively.

4.1.

Temporal Evolution of the Light Propagation

In the forward problem we use Eq. 4 to model the light propagation and resulting fluorescence TPSF given a priori optical properties of a turbid medium and the fluorophore inclusion properties. Figure 2 shows the evolution of the fluorescence TPSF at each propagation step with normalized maximum intensity to highlight the increasing temporal convolution form. To validate the forward model, we used the optical properties from the phantom experiment with a Cy5 fluorophore inclusion (i.e., ${\mu }_a=3.0\times {10}^{-4}{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.3{\mathrm{mm}}^{-1}$ , $\tau =1.5\mathrm{ns}$ , and $d=1\mathrm{mm}$ ). Rather than attempting to calibrate the experiment intensity units to model intensity units, we simply display the normalized maximum intensity experimental data in Fig. 2 and demonstrate that the forward model can predict the correct shape of the experimental fluorescence TPSF. Moreover, since the model assumes a point fluorophore and experimentally we have a finite fluorescent inclusion, it is not strictly possible to model the absolute fluorophore concentration, although we can still compare fluorophore concentration in arbitrary units within this point-model approximation. From Fig. 2 it can also be seen that the time decay of the fluorescence TPSF is largely dominated by the lifetime decay of the fluorescent inclusion. However, the so-called effective lifetime ${\tau }_{\mathrm{eff}}$ indicating the time decay of a fluorescence TPSF is not the same as the actual lifetime $\tau$ of a fluorescent inclusion due to the convolution of its lifetime decay with two Green’s functions and the IRF. For more detailed analysis, the increase in the decay of the fluorescence TPSF is evaluated by fitting it with a monoexponential function from 80% peak intensity to 20% peak intensity as the photons propagate from source to detector. The first Green’s function $\left[f_1\left(t\right)\right]$ , i.e., the propagation of light from source $\left({\mathbf{r}}_s\right)$ to a fluorescent inclusion $\left(\mathbf{r}\right)$ in a turbid medium, shows a rapid decay of only $33\mathrm{ps}$ , which does not contribute significantly to the effective lifetime of the resultant fluorescence TPSF, which is $1.98\mathrm{ns}$ . After photons excite the fluorescent inclusion to cause fluorescence emission, the decay of the consequent function $[f_1\left(t\right)\ast f_2\left(t\right)]$ is $1.64\mathrm{ns}$ , which is a significant component of the effective lifetime. When convoluted with another Green’s function $\left[f_3\left(t\right)\right]$ , there is a slight increase in the decay of the consequent function $[f_1\left(t\right)\ast f_2\left(t\right)\ast f_3\left(t\right)]$ to $1.78\mathrm{ns}$ . Finally, convolution with the system IRF $\left[f_{\mathrm{IRF}}\right]$ adds further increase, resulting in the fluorescence TPSF $[f_1\left(t\right)\ast f_2\left(t\right)\ast f_3\left(t\right)\ast f_{\mathrm{IRF}}\left(t\right)]$ , whose decay is $1.98\mathrm{ns}$ and is the effective lifetime.

Fig. 2

Temporal evolution of the maximum intensity normalized fluorescence TPSF at each step during photon migration from source to detector.

4.2.

Effect of Background Optical Properties on the Fluorescence Temporal Point-Spread Function

We assumed the background optical properties were known accurately a priori. However, in vivo we have to estimate the optical properties of the mouse where inaccuracy can introduce error in recovering $n$ , $\tau$ , and $d$ of the fluorophore. To mimic this condition, we used the forward model to generate fluorescence TPSFs for a variety of background optical properties (5% change in ${\mu }_a$ or ${\mu }_s^{\prime }$ from initial values). We then fitted these fluorescence TPSFs with the model fixing the background optical properties to the initial values and investigated the impact on the recovered fitted values of $n$ , $\tau$ , and $d$ . We found that there was negligible error in $\tau$ and $d$ (0.5%). Interestingly, Fig. 3 shows the fluorescence TPSFs normalized to maximum intensity, where it is observed that the temporal shape of the fluorescence TPSF for the different values of ${\mu }_a$ [Fig. 3a] and ${\mu }_s^{\prime }$ [Fig. 3b] is similar. Since previous work⁸ has shown that $\tau$ and $d$ are strongly correlated with the temporal shape of the fluorescence TPSF, i.e., ${\tau }_{\mathrm{eff}}$ and $t_{\max }$ , it is expected that similar $\tau$ and $d$ values were recovered. However, varying the background optical properties also influences the absolute area of the fluorescence TPSF, which was found to induce a 10% error in $n$ that essentially scales the intensity signal without changing the temporal shape. Nevertheless, recovery of fluorophore concentration within 10% would still be a useful result, given that the background optical properties can be estimated a priori within 5%.

Fig. 3

Maximum intensity normalized fluorescence TPSFs for varying background optical properties for (a) ${\mu }_a=3.0\times {10}^{-5}–3.0\times {10}^{-3}{\mathrm{mm}}^{-1}$ and (b) ${\mu }_s^{\prime }=1.2\text{to}1.4{\mathrm{mm}}^{-1}$ .

4.3.

Phantom Analysis

In this section we analyze experimental phantom data from Cy5 and Cy7 fluorescent inclusions submerged in a turbid medium. The background optical properties of the phantom were as follows: ${\mu }_a=3.0\times {10}^{-4}{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.3{\mathrm{mm}}^{-1}$ at $\lambda =635\mathrm{nm}$ (for Cy5), and ${\mu }_a=2.0\times {10}^{-3}{\mathrm{mm}}^{-1}$ , and ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$ at $\lambda =760\mathrm{nm}$ (for Cy7).¹⁵ The eXplore-Optix-MX2 scanned the phantom and acquired fluorescence TPSFs as described before. Model fitting involved first averaging the fluorescence TPSFs from a region of interest covering the fluorescent pellet to extract an average value of $n$ , $\tau$ , and $d$ with a priori background optical properties. Using the average values as an initial guess, we then fitted the model [Eq. 5] on a pixel-by-pixel basis. The difference between the measured and calculated fluorescence TPSF in each pixel was minimized using a least-squares nonlinear optimization routine, i.e., the Levenberg-Marquardt algorithm.^{16, 17} For example, Fig. 4 shows a final fit result for a given pixel from the Cy5 fluorescent inclusion after 88 iterations from initial average parameters with a priori knowledge of the background optical properties ( ${\mu }_a=3.0\times {10}^{-4}{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.3{\mathrm{mm}}^{-1}$ ), as well as $t_{\mathrm{IRF}}=1.00\mathrm{ns}$ , and ${\sigma }_{\mathrm{IRF}}=0.25\mathrm{ns}$ at $\lambda =635\mathrm{nm}$ . Recovered fitted values of $n=6.5\times {10}^4{\mathrm{mm}}^{-3}$ , $\tau =1.5\mathrm{ns}$ , and $d=1.3\mathrm{mm}$ were obtained.

Fig. 4

Example fit result from Cy5 inclusion using the algorithm described in Eq. 5 with a priori values of ${\mu }_a=3.0\times {10}^{-4}{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.31{\mathrm{mm}}^{-1}$ , $t_{\mathrm{IRF}}=1.00\mathrm{ns}$ , and ${\sigma }_{\mathrm{IRF}}=0.25\mathrm{ns}$ at $\lambda =635\mathrm{nm}$ . Recovered fit values are evaluated as $n=6.5\times {10}^4{\mathrm{mm}}^{-3}$ , $\tau =1.52\mathrm{ns}$ , and $d=1.3\mathrm{mm}$ .

Displaying the recovered fitted values of $n$ , $\tau$ , and $d$ for each scanned pixel generated corresponding 2-D images of these parameters. Since the model assumes a point fluorophore at depth $d$ directly below the surface, we recover appropriate values for the inclusion as we scan across it. However, as we scan away from the inclusion, the model breaks down and we no longer recover appropriate values. In essence, the model corrects for optical propagation in the $z$ axis but not in the $x\text{-}y$ scanning plane where optical scattering increases the apparent size of the inclusion. Figure 5 show images of $n$ (top row), and $d$ (bottom row) of the Cy5 fluorescent inclusion submerged at $1\mathrm{mm}$ (left) and $5\mathrm{mm}$ (right), respectively. Figure 6 shows the corresponding images for the Cy 7 fluorescent inclusion. In Fig. 5 (top), it is seen in the vicinity of the inclusion that similar values of $n(\sim 1.0\times {10}^6)$ are recovered at 1- and $5\text{-}\mathrm{mm}$ submersion as expected, since the fluorophore concentration remained constant, even though the cw intensity diminished with increasing depth. In Fig. 5 (middle), it is seen in the vicinity of the inclusion that similar values of $\tau (\sim 1.5\mathrm{ns})$ are recovered at 1- and $5\text{-}\mathrm{mm}$ submersion as expected, since the fluorophore lifetime remained constant even though the effective lifetime increased with increasing depth. In Fig. 5 (bottom), it is seen in the vicinity of the inclusion that different values of $d$ ( $\sim 1.3$ and $5.4\mathrm{mm}$ ) are recovered at 1- and $5\text{-}\mathrm{mm}$ submersion as expected, since the fluorophore depth increased. Note the submersion depth is defined as the perpendicular distance to the top surface of the fluorescent inclusion. Given that we have a point model and a finite inclusion, it is expected to slightly overestimate $d$ , as discussed elsewhere.⁸ Similar results were found for the Cy 7 fluorescent inclusion as seen in Fig. 6, where fitted values of $n(\sim 1.2\times {10}^6)$ , $\tau (\sim 1\mathrm{ns})$ , and $d$ ( $\sim 1.3$ and $5.2\mathrm{mm}$ ) were recovered.

Fig. 5

Images of $n$ (top row), $\tau$ (middle row), and $d$ (bottom row) for shallowly (left) and deeply (right) embedded Cy5 inclusion.

Fig. 6

Images of $n$ (top row), $\tau$ (middle row), and $d$ (bottom row) for shallowly (left) and deeply (right) embedded Cy7 inclusion.