1.

Introduction

Fluorescence guidance in neurosurgery using 5-aminolevulinic acid (ALA)-induced protoporphyrin IX (PpIX) is gaining popularity as a surgical adjunct to improve the extent of tumor resection.^1–3 A randomized controlled phase III clinical trial by Stummer and colleagues has shown greater extent of resection and significantly higher 6-month progression-free survival in patients with high-grade glioma (GBM) undergoing fluorescence guidance compared to conventional white-light-guided surgery.^4,5 Methods to quantify PpIX concentration in tissue using a fiber-optic “point” probe and a spectrally constrained model of light propagation have been shown to improve tumor detection with enhanced accuracy (87%) compared to the subjective visual fluorescence assessment (66%) used in the Stummer studies.^6,7 To better integrate this concept into surgical guidance, this work explores the development of a model-based wide-field fluorescence imaging method integrated into a surgical microscope. Our goal is to develop quantitative, depth-resolved imaging for fluorescence guidance during neurosurgery.

To date, ALA-induced PpIX fluorescence guidance has been based on violet/blue light to excite the Soret-band of the molecule (around $\lambda \sim 405\mathrm{nm}$), which limits the penetration of the probing signal to a few hundred microns at best, potentially leaving residual tumor undetected when located below the surgical surface. Fortunately, the excitation spectrum of PpIX also possesses several longer wavelength excitation bands up to around $\lambda \sim 633\mathrm{nm}$ (Q-bands) and a broad emission spectrum extending up to $\sim 720\mathrm{nm}$ that could be used to stimulate (and detect) emissions from subsurface tumor. Despite the potential benefits during surgery, efforts to detect fluorescence at depth using model-based techniques are still in an early stage of development with promising preliminary results reported by Kim et al.⁸ and Leblond et al.⁹

Spatial frequency domain imaging (SFDI) is a relatively new method that uses red or near-infrared light together with CCD detection to acquire wide-field images of tissue. Sinusoidal patterns of varying spatial frequencies are projected onto tissue, and the attenuation in the amplitude of the spatial patterns at each CCD pixel is fitted to a model of light propagation in order to form images of the tissue optical properties.^10,11 SFDI has recently been used to quantify PpIX concentrations in vivo by accounting for tissue optical properties.¹² Spatially modulated excitation also encodes the depth from which measured optical signals originate because high spatial frequency patterns penetrate less deeply into tissue than their lower frequency counterparts because of multiple light scattering. This frequency encoding enables depth information to be recovered through three-dimensional (3-D) optical tomography.¹³

Most SFDI research has focused on generating two-dimensional (2-D) images by fitting the reflectance measure at each pixel of a CCD array to a homogeneous model of light propagation. With the approach, SFDI has been applied to quantitatively image stroke,¹⁴ brain injury,¹⁵ cortical spreading depression,¹⁶ and layered structures in skin.¹⁷ SFDI techniques have also been used to localize subsurface inhomogeneities,¹⁸ and to determine tissue optical absorption and scattering properties over a broad spectral range.¹⁹ Recently, several groups have reconstructed 3-D tomographic images using sinusoidal illumination patterns,^13,20,21 but these studies were limited to measuring the optical properties of tissue simulating phantoms and did not include in vivo or fluorescence measurements. While we have performed fluorescence SFDI, these preliminary studies did not involve tomographic image reconstruction and were limited to phantoms and/or skin.^12,22 Likewise, while fluorescence optical tomography with wide-field illumination has been reported for phantom studies^23,24 and whole-body mouse imaging,²⁵ it has not been reported either for spatial frequency domain illumination or brain imaging.

In this study we implement, for the first time to our knowledge, fluorescence optical tomography using spatial frequency domain illumination. We then demonstrate the ability of SFDI to image GBM tumors based on endogenous contrast from hemoglobin and the tissue-scattering parameters as well as exogenous contrast from ALA-induced PpIX fluorescence. We present experiments in tissue simulating phantoms with embedded GBM spheroids, and in an in vivo GBM mouse model. We use 2-D SFDI to image endogenous optical properties and observe that the reduced scattering coefficient (${\mu }_s^{\prime }$) exhibits significant contrast between normal brain and tumor. We then successfully implement 3-D SFDI fluorescence tomography to localize tumors through ALA-induced PpIX fluorescence.

2.

Methods

2.2.

Data Acquisition

A schematic of the instrument appears in Fig. 1. We imaged reflectance and fluorescence by illuminating a liquid crystal on silicon (LCOS) device (HEO 1080p, HOLOEYE Photonics AG) with a 20-mW, $\sim 633\text{-}\mathrm{nm}$ HeNe laser (1145AP, JDS Uniphase). The LCOS projected images on the sample consisted of sinusoidal intensity patterns of the form

Eq. (1)

$$S_i(\mathbf{\rho })=\frac{1}{2}[1+\cos (\mathbf{k}·\mathbf{\rho }+{\varphi }_i)],$$

where $\mathbf{\rho }$ represents a location on the sample surface, $\mathbf{k}$ is the spatial projection wave number, ${\varphi }_i$ is the phase of the spatial projection. The spatial projection frequency $f=2\pi |\mathbf{k}|$ was varied from $f=0{\mathrm{mm}}^{-1}$ (no spatial modulation) to $f=0.3{\mathrm{mm}}^{-1}$ in intervals of $f=0.05{\mathrm{mm}}^{-1}$. At each spatial frequency, the intensity pattern was serially projected at three different phases (${\varphi }_i=0$, 120, and 240 deg), and the remitted light resulting from each projection was detected by a CCD array (Cascade 512B, Roper Scientific). The average excitation intensity at the tissue surface for $f=0{\mathrm{mm}}^{-1}$ was $\sim 1\mathrm{mW}/{\mathrm{cm}}^2$. In fluorescence mode, the excitation light was blocked at the detector with a 650-nm longpass filter ($\text{peak}\text{transmission}=90\%$, out-of-band $\mathrm{OD}>5$) followed by a 680-nm bandpass filter [full width at half maximum $(\mathrm{FWHM})=20\mathrm{nm}$, $\text{peak}\text{transmission}=70\%$, out-of-band $\mathrm{OD}>3$]. To image endogenous contrast, light from a quartz tungsten halogen lamp (CF1000, Illumination Technologies) was filtered serially at four wavelengths (680, 720, 800, and 820 nm) and coupled to the LCOS via a 5-mm-core-diameter liquid light guide. Orthogonal linear polarizers were incorporated to eliminate glare from specular reflections and were used in both the endogenous and fluorescence scans. Data acquisition took $\sim 30\mathrm{s}$ (endogenous contrast), and $\sim 100\mathrm{s}$ (fluorescence contrast). We expect to be able to reduce these times by use of brighter sources and optimized optics.

Fig. 1

Schematic of the SFDI instrument. Filters were removed when imaging endogenous contrast.

2.3.

Data Analysis

The three images ($I_1$, $I_2$, $I_3$) at the different phases acquired for each projected spatial frequency were combined according to Ref. 13:

Eq. (2)

$$I_{\mathrm{tot}}(\mathbf{\rho })=\frac{1}{3}[2I_2(\mathbf{\rho })-I_1(\mathbf{\rho })-I_3(\mathbf{\rho })]+\frac{i}{\sqrt{3}}[I_1(\mathbf{\rho })-I_3(\mathbf{\rho })].$$

As can be seen by substituting the source term described by Eq. (1) into Eq. (2), $I_{\mathrm{tot}}(\rho )$ represents the reflectance expected from a source of the form $S_{\mathrm{tot}}(\mathbf{\rho })\sim \exp (i\mathbf{k}·\mathbf{\rho })$. To produce 2-D maps of the absorption coefficient, ${\mu }_a$, and the reduced scattering coefficient, ${\mu }_s^{\prime }$, at each wavelength, the amplitude of $I_{\mathrm{tot}}(\rho )$ at each CCD pixel was fitted to a semi-infinite, homogeneous, radiative transport model of light propagation aided by a rapid look-up table approach.¹⁰ Tissue concentrations of oxyhemoglobin (${\mathrm{ctHbO}}_2$) and deoxyhemoglobin (ctHb) were determined through spectral fitting of the absorption values using the known extinction coefficients for the two chromophores (assuming no other significant absorbers were present).

We reconstructed 3-D images of the concentration of PpIX by adapting our previously published method for tomographic imaging of absorption using SFDI.¹³ For continuous wave illumination, the measured emission light at position ${\rho }_d$ due to illumination wavenumber $\mathbf{k}$ is related to the fluorophore concentration by

Eq. (3)

$$I_{\mathrm{tot}}(\mathbf{k},{\mathbf{\rho }}_d)=\int {\mathrm{d}}^{3}r{I}_{\mathrm{ex}}(\mathbf{k},\mathbf{r})G_{\mathrm{em}}(\mathbf{r},{\mathbf{r}}_d)\epsilon ({\lambda }_{\mathrm{ex}})\eta ({\lambda }_{\mathrm{ex}},{\lambda }_{\mathrm{em}})C(\mathbf{r}),$$

where $I_{\mathrm{ex}}(\mathbf{k},\mathbf{r})$ is the photon density of excitation light in the medium, $\epsilon ({\lambda }_{\mathrm{ex}})$ is the extinction coefficient of the fluorophore at the excitation wavelength, $\eta ({\lambda }_{\mathrm{ex}},{\lambda }_{\mathrm{em}})$ is the fluorescence quantum yield ($<1\%$ for PpIX), and $C(\mathbf{r})$ is the fluorophore concentration. The Green’s function $G_{\mathrm{em}}(\mathbf{r},{\mathbf{r}}_d)$ obeys the diffusion equation at the emission wavelength

Eq. (4)

$$[-\nabla ·D\nabla +c{\mu }_a]G(\mathbf{r},{\mathbf{r}}^{\prime })=\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

subject to the boundary condition

Eq. (5)

$$G(\mathbf{r},{\mathbf{r}}^{\prime })+\ell \widehat{\mathbf{n}}·\nabla G(\mathbf{r},{\mathbf{r}}^{\prime })=0.$$

Here $D=3/c\ell *$ is the diffusion coefficient, $\ell$ is the extrapolation length, and $\widehat{\mathbf{n}}$ is the outward-pointing unit normal. For the semi-infinite geometry, the Green’s function can be decomposed into plane waves²⁸ according to

Eq. (6)

$$G(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{{(2\pi )}^2}\int {\mathrm{d}}^{2}qg(\mathbf{q},z,z^{\prime })\exp \{i[\mathbf{q}·({\mathbf{\rho }}^{\prime }-\mathbf{\rho })]\},$$

where

Eq. (7)

$$g(\mathbf{q},z,z^{\prime })=\frac{\ell }{D}\frac{\exp [-Q(\mathbf{q})|z-z^{\prime }|]}{Q(\mathbf{q})\ell +1},$$

when either $z$ or $z\prime$ is located on the surface, and

Eq. (8)

$$Q(\mathbf{q})={(c{\mu }_a/D+{\mathbf{q}}^2)}^{1/2}.$$

Equations (6) and (7) allow us to represent the photon density of the exitation light in the tissue appearing in Eq. (3) as

Eq. (9)

$$I_{\mathrm{ex}}(\mathbf{k},\mathbf{r})=g_{\mathrm{ex}}(\mathbf{k},0,z)\exp [i\mathbf{k}·\mathbf{\rho }].$$

Substituting Eqs. (6), (7), and (9) into Eq. (3) and taking a Fourier transform with respect to the detector positions on the surface result in

Eq. (10)

$$I_{\mathrm{tot}}(\mathbf{k},{\mathbf{q}}_d)=\epsilon ({\lambda }_{\mathrm{ex}})\eta ({\lambda }_{\mathrm{ex}},{\lambda }_{\mathrm{em}})\int \mathrm{d}z\kappa (\mathbf{k},{\mathbf{q}}_d,z)\int {\mathrm{d}}^2\rho \exp [i\mathbf{q}·\mathbf{\rho }]C(\mathbf{r}),$$

where

Eq. (11)

$$\kappa (\mathbf{k},{\mathbf{q}}_d,z)=\frac{{\ell }_{\mathrm{ex}}{\ell }_{\mathrm{em}}}{D_{\mathrm{ex}}{D}_{\mathrm{em}}}\frac{\exp \{-[Q_{\mathrm{ex}}(\mathbf{k})+Q_{\mathrm{em}}({\mathbf{q}}_d)]z\}}{[Q_{\mathrm{ex}}(\mathbf{k}){\ell }_{\mathrm{ex}}+1][Q_{\mathrm{em}}({\mathbf{q}}_d){\ell }_{\mathrm{em}}+1]},$$

and

Eq. (12)

$$I_{\mathrm{tot}}(\mathbf{k},{\mathbf{q}}_d)=\int {\mathrm{d}}^2\rho \exp (i{\mathbf{q}}_d·\mathbf{\rho })I_{\mathrm{tot}}(\mathbf{k},{\mathbf{\rho }}_d).$$

The inverse of Eq. (10) is used to calculate the fluorophore distribution. It appears as

Eq. (13)

$$C(\mathbf{r})=\frac{1}{\epsilon ({\lambda }_{\mathrm{ex}})\eta ({\lambda }_{\mathrm{ex}},{\lambda }_{\mathrm{em}})}\int \frac{d^{2}q}{{(2\pi )}^2}\exp (-i\mathbf{q}·\mathbf{\rho })\sum _{m,n}{\kappa }^{*}({\mathbf{k}}_m,{\mathbf{q}}_d,z)M_{m,n}^{-1}(\mathbf{q})I_{\mathrm{tot}}({\mathbf{k}}_n,{\mathbf{q}}_d),$$

where

Eq. (14)

$$M_{m,n}(\mathbf{q})={\int }_0^{\infty }dz\kappa ({\mathbf{k}}_m,{\mathbf{q}}_d,z){\kappa }^{*}({\mathbf{k}}_n,{\mathbf{q}}_d,z).$$

In these formulas, the Fourier variable $\mathbf{q}$ is related to the projection wave number $\mathbf{k}$ and the Fourier variable conjugate to the detector locations ${\mathbf{q}}_d$ by $\mathbf{q}=\mathbf{k}+{\mathbf{q}}_d$. The subscripts ex and em indicate whether the optical properties (${\mu }_a$ and ${\mu }_s^{\prime }$) for the respective variable refer to the excitation or emission wavelength ($\sim 633\mathrm{nm}$ and $680\mathrm{nm}\pm 20\mathrm{nm}$, respectively). The average optical properties at 633 nm (excitation) and 680 nm (emission) obtained from the endogenous SFDI were used as inputs in Eq. (5) for the tomographic fluorescence reconstructions. The inverses of the M matrices [see Eq. (14)] are needed in Eq. (13). We take the inverse of each of these matrices by solving its eigenvalue problem and setting all eigenvalues below a threshold to zero. For this study the threshold was kept constant at ${10}^{-4}$. We solved for the fluorophore concentration over a $50\times 50$ grid of points in the Fourier domain with $q_{\max }=2.5\pi {\mathrm{mm}}^{-1}$ and $\delta q=0.1\pi {\mathrm{mm}}^{-1}$. Image reconstruction required only $\sim 10\mathrm{s}$ on a single 2.80-GHz processor. A flowchart of the tomographic reconstruction algorithm appears in Fig. 2.

Fig. 2

Flowchart of the three-dimensional fluorescence optical tomography algorithm.

3.

Results

As an initial test of the ability to image relative PpIX concentrations in turbid media, two GBM spheroid phantoms were evaluated. In the first phantom the spheroid was located with its top surface level with the surface of the phantom, whereas in the second phantom the spheroid was embedded such that its top surface was located $\sim 1.5\mathrm{mm}$ beneath the phantom. After imaging, the phantom containing the subsurface spheroid was cut in half to confirm its location [Fig. 3(a)]. In the reconstructed images, the maximum PpIX concentration of the superficial spheroid [Fig. 3(b)] was reconstructed at a depth of 1 mm, while the maximum for the subsurface spheroid [Fig. 3(c)] was at a depth of 2 mm. Although both spheroids were approximately the same size ($\sim 1.5\mathrm{mm}$), the image of the superficial spheroid had a lateral FWHM of 2.4 mm compared to 4.7 mm for the deeper spheroid. The volume of contrast enhancement in the reconstructed image with concentration greater than 50% of the maximum was also less for the superficial ($13.7{\mathrm{mm}}^3$) than for the deeper ($66.8{\mathrm{mm}}^3$) spheroid. However, the reconstructed maximum concentration of the superficial spheroid was 3.6 times that of the deeper spheroid, such that the estimated total amount of PpIX in the two images agreed within 30%.

Fig. 3

Fluorescence tomography reconstructions of PpIX concentrations for the GBM spheroid containing phantoms. (a) Cross-sectional image of the subsurface spheroid taken after the phantom was cut in half. (b) Reconstructed image of the phantom with the spheroid at the surface. (c) As in (b) for the subsurface spheroid. Each image in (b and c) corresponds to the reconstructed PpIX concentration at the depth indicated below the frame.

In the mouse images of the reduced scattering coefficient and hemoglobin concentration through the intact skull [Fig. 4(a)] and with the brain exposed [Fig. 4(b)], scattering contrast between tumor and normal tissue (taken as an identically sized region adjacent to the tumor) was $>1.5$ at all wavelengths [Fig. 4(c)]. Fitting scattering values to a power law of the form ${\mu }_s^{\prime }=A{(\lambda /{\lambda }_0)}^{-b}$ with ${\lambda }_0=750\mathrm{nm}$ resulted in A and b values that were greater for the tumor. As shown in Fig. 4(d), deoxyhemoglobin concentration, ctHb, was significantly higher in tumor than in normal tissue, likely because of metabolic differences such as angiogenesis and hypoxia. Oxyhemoglobin concentration values, ${\mathrm{ctHbO}}_2$, were similar for tumor and normal brain tissue. These results are summarized in Table 1.

Fig. 4

Two-dimensional images of the reduced scattering coefficient, oxyhemoglobin concentration, and deoxyhemoglobin concentration in the mouse with (a) intact skull and (b) exposed brain. Plots of the scattering coefficient, ${\mu }_s^{\prime }$, (c) and (d) concentrations of Hb and ${\mathrm{HbO}}_2$ for tumor and normal regions are also shown. Error bars denote the standard deviation of pixel values in the region of interest.

Table 1

Measured in vivo values in the mouse brain for tumor and normal tissue (Shown are the means and standard deviations of pixel values in the tumor and normal regions of interest).

Figure 5(a) shows a picture of the skull illuminated with spatially modulated 633-nm light ($f=0.25{\mathrm{mm}}^{-1}$). The dotted red line denotes the region of interest used in the fluorescence image reconstruction. In the reconstructed images of PpIX concentration, the tumor appears slightly deeper with an intact skull [Fig. 5(b)] than with the skull removed [Fig. 5(c)], with the maximum intensity occurring 1.5 mm and 1.0 mm below the surface, respectively (actual skull thickness $\sim 0.5\mathrm{mm}$). The presence of the skull decreased the maximum concentration in the image by 8%, increased the FWHM from 1.1 to 2.1 mm, and increased the apparent tumor volume from 4.5 to $5.5{\mathrm{mm}}^3$. Similar to the phantom experiments, the effects of increased size and decreased concentration offset each other such that the total estimated amount of PpIX remained constant to within 10%.

Fig. 5

(a) Photograph of the skull in which the illumination pattern is visible. The dotted red line denotes the region of interest for the 3-D reconstructions. (b) Reconstructed images of PpIX concentration with the skull intact. (c) As in (b) with the brain surgically exposed. The corresponding reconstructed image depth is indicated below each frame.

4.

Discussion

One feature common to all of the reconstructions in this study is that deeper objects appear larger and with less contrast than superficial ones. This is the result of a partial volume effect that occurs because the interrogation volume of each measurement is larger than the features we are trying to reconstruct. As a result, the reconstructed features have values somewhere between those of their actual values and the background. The effect is more pronounced for deeper structures because the sensitivity functions are peaked near the surface and become broader with depth. We have analyzed the partial volume effect for SFDI in a previous publication where we showed that tomographic image reconstruction ameliorated (but did not eliminate) this effect in comparison to our 2-D imaging methods.¹³

We do not expect this partial volume effect to prevent the use of SFDI tomography for neurosurgical resection of tumors. Fluorescence-assisted surgical resection is guided by the qualitative presence or absence of fluorescent signal. Thus, the exact peak concentration is of less importance than a measurement of the fluorophore depth or the total amount of fluorophore present. Knowledge of fluorophore depth could improve the efficiency of surgical resection by allowing a surgeon to more confidently dissect down to tumor below the cortical surface; then proceed with more aggressive tumor resection while far from the deep margin; and finally convert to slower, more meticulous removal when the deep interface is near. The increase in size of the reconstructions with depth is expected to interfere with accurate estimates of deep tumor extent. However, as tumor can only be removed from the surface of the operative field, the higher volume estimates at depth are of less concern because they are expected to decrease as tumor is removed and overall depth is reduced.

It may be possible to obtain more accurate reconstructions by using algorithms that allow for inhomogeneous tissue optical properties. Such algorithms could be implemented with finite-element or Monte Carlo methods to model light transport. However, one of the main advantages of our reconstruction method is its speed. This speed is necessary to make surgical applications feasible. The current reconstructions use $\sim {10}^6$ data points and take only $\sim 10\mathrm{s}$. In addition, they could be made faster by parallel computing, by pre-calculating the eigenvalues/functions of the $M$ matrices, or by reconstructing over a reduced bandwidth. We are able to attain this speed because we use homogeneous tissue optical properties and a planar geometry to model light transport. These assumptions make the Green’s function for the diffusion equation translationally invariant and enable us to write Eq. (3) in block diagonal form [see Eq. (10)] such that it can be solved rapidly. Thus while we plan to explore the potential benefits of modeling light transport through inhomogeneous tissue, we expect that for applications requiring rapid feedback (e.g., surgery), the algorithm proposed here will be better suited.

5.

Conclusion

In summary, this proof-of-principle study demonstrates the ability of multiwavelength, multifrequency SFDI tomography to localize small ($\sim \mathrm{mm}$) subsurface ($\sim 2$ to 3 mm deep) tumors in tissue-simulating phantoms and in vivo using ALA-induced PpIX fluorescence. These results suggest that SFDI may offer a wide-field fluorescence imaging approach that accounts for intrinsic tissue optical properties, characterizes subsurface nests of tumor cells, and potentially provides contrast for improved intraoperative neurosurgical guidance. Future studies will provide more detailed insight into factors that limit SFDI accuracy and performance, and explore the combination of image guidance with probe-based methods.^6,8,9