## Abstract

A promising method to incorporate tissue structural information into the reconstruction of diffusion-based fluorescence imaging is introduced. The method regularizes the inversion problem with a Laplacian-type matrix, which inherently smoothes pre-defined tissue, but allows discontinuities between adjacent regions. The technique is most appropriately used when fluorescence tomography is combined with structural imaging systems. Phantom and simulation studies were used to illustrate significant improvements in quantitative imaging and linearity of response with the new algorithm. Images of an inclusion containing the fluorophore Lutetium Texaphyrin (Lutex) embedded in a cylindrical phantom are more accurate than in situations where no structural information is available, and edge artifacts which are normally prevalent were almost entirely suppressed. Most importantly, spatial priors provided a higher degree of sensitivity and accuracy to fluorophore concentration, though both techniques suffer from image bias caused by excitation signal leakage. The use of spatial priors becomes essential for accurate recovery of fluorophore distributions in complex tissue volumes. Simulation studies revealed an inability of the “no-priors” imaging algorithm to recover Lutex fluorescence yield in domains derived from T1 weighted images of a human breast. The same domains were reconstructed accurately to within 75% of the true values using prior knowledge of the internal tissue structure. This algorithmic approach will be implemented in an MR-coupled fluorescence spectroscopic tomography system, using the MR images for the structural template and the fluorescence data for region quantification.

© 2007 Optical Society of America

## 1. Introduction

Imaging the spatial distribution of fluorescence activity at depth in tissue is a challenging problem that has yet to impact clinical practice. The most popular approach is diffuse optical fluorescence tomography (DOFT), a model-based method which approximates photon propagation as a diffuse field and computationally matches model parameters to measured boundary data. The theoretical framework is well established [1–5] and feasibility studies have demonstrated fluorescence yield imaging in various phantom geometries [6–9]. Much of the most recent research has focused on high-throughput, full body imaging of small animals [10, 11] but no human images have been published to date. DOFT produces low resolution images compared to standard clinical modalities such as x-ray CT and MRI due to the highly scattered photon fields and relatively sparse measurement sampling of the tissue volume. The reconstruction problem is ill-posed and underdetermined and large heterogeneity in tissue optical properties caused by complex tissue morphology further challenges the imaging algorithm. Sensitivity drops with increasing depth resulting in a non-linear responsivity across the imaging field, making it more difficult to image larger tissue volumes. Experiments in larger volumes that more closely resemble a human breast typically consider unrealistically high fluorophore contrasts and simple geometries [9, 12]. Though these studies are important for advancing the understanding of the modality, improving image resolution and contrast sensitivity is critical for identifying a clinical role for DOFT.

DOFT is an extension of the more widely studied diffuse optical tomography (DOT) which suffers similarly from low resolution and depth-dependent contrast sensitivity in the absence of spatial or spectral prior information to guide the solution. However, methods to incorporate highly-resolved anatomical data obtained from standard clinical modalities have improved the quantification of the optically derived images [13–15]. These hybrid approaches lead to a conceptually new application of optical tomography, one in which the highly resolved imaging system provides a structural template upon which volumetric optical spectroscopic images are constructed. This framework may be applied to either absorption and scatter spectroscopy, or fluorescence spectroscopy. Additional challenges specific to the fluorescence case include lower signal intensity, excitation source contamination of the fluorescence emission measurements due to filter leakage, and tissue optical property effects on both the excitation and emission photon propagation. Applying spatial guidance techniques to the more complicated DOFT problem may yield even larger gains in imaging capability.

This paper introduces a method to dramatically improve fluorescence imaging at depth in tissue by coupling MRI and fluorescence tomography. Tissue structural information determined from standard T1 and T2 MR images is encoded as a spatial filter in the DOFT reconstruction algorithm and used to guide the recovery of fluorescence activity. In this implementation, reconstruction parameters are loosely grouped into regions based on tissue-type determined from the MR images but are permitted to update independently, giving rise to the term “soft” spatial priors [16]. The algorithm is tested with phantom data of Lutetium Texaphyrin, a near-infrared photosensitizer shown to accumulate preferential in a variety of malignancies [17–20], recorded from a newly developed MR-coupled spectroscopy-based fluorescence tomography system. A simulation study based on realistic geometries generated from MR images of a normal human breast serves as an initial illustration of expected improvements provided by the spatial priors approach in complex domains.

## 2. Theory

Assuming highly diffuse media, the transport of light at a given modulation frequency ω in the presence of fluorescence generated by an external source at the excitation wavelength (λ_{x}) of the fluorescing agent, is modeled by two diffusion equations, where the solution of the first equation provides the driving source term of the second [5]:

where subscripts *x* and *m* represent the excitation and emission fluence at wavelengths λ_{x} and λ_{m}, respectively. The intrinsic optical parameters *μ _{ax,m}* and ${\mu}_{{s}_{x,m}}^{\prime}$ are the absorption and reduced scattering coefficients respectively,

*q*

_{0}(

*r*,

*ω*) is an isotropic source and Φ

_{x,m}(

*r*,

*ω*) is the photon fluence rate at position

*r*. The diffusion coefficient is given by ${\kappa}_{x,m}=\frac{1}{3\left({\mu}_{{a}_{x,m}}+{\mu}_{{s}_{x,m}}^{\prime}\right)}$ and

*c(r)*is the speed of light in the medium at any point, defined by

*c*, where

_{o}/n(r)*n(r)*is the index of refraction at the same location and

*c*is the speed of light in vacuum. The fluorescence parameters are the lifetime

_{o}*τ*(

*r*) and the fluorescence yield

*ημ*(

_{af}*r*), the latter a product of the fluorophore’s quantum efficiency

*η*and its absorption coefficient

*μ*(

_{af}*r*).

The most appropriate description of the air-tissue boundary is derived with an index-mismatched type III condition (also known as Robin or mixed), in which some fraction of the fluence at the external boundary of the tissue exits and does not return. The flux leaving the external boundary is equal to the fluence rate at the boundary weighted by a factor that accounts for the internal reflection of light back into the tissue. This relationship is described in the following equation:

where *ξ* is a point on the external boundary, and *A* depends upon the relative refractive index (RI) mismatch between tissue Ω and air. *A* can be derived from Fresnel’s law:

where *θ _{C}* = arcsin(

*n*/

_{AIR}*n*

_{1}), the angle at which total internal reflection occurs for photons moving from region Ω with RI

*n*

_{1}to air with RI

*n*

_{AIR}, and ${R}_{0}=\frac{{\left(\frac{{n}_{1}}{{n}_{AIR}}-1\right)}^{2}}{{\left(\frac{{n}_{1}}{{n}_{AIR}}+1\right)}^{2}}$. At the external boundaries,

*n*= 1.

_{AIR}#### 2.1 Finite element implementation:

When the refractive index is homogenous (assumed to be 1.33 [21]), the finite element discretization of a volume Ω can be obtained by subdividing the domain into D elements joined at V vertex nodes. In the finite element method (FEM), Φ_{x,m}(**r**) is approximated by the piecewise continuous polynomial function ${\Phi}_{x,m}^{h}\left(r,\omega \right)={\sum}_{i}^{V}{\Phi}_{x,{m}_{i}}{u}_{i}\left(r\right){\Omega}^{h}$, where Ω^{h} is a finite dimensional subspace spanned by basis functions {*u _{i}*(

*r*);

*i*= 1⋯

*V*} chosen to have limited support. The problem of solving for Φ

^{h}

_{x,m}becomes one of sparse matrix inversion, and in this work, a bi-conjugate gradient stabilized solver is used. As developed previously [22, 23], the coupled diffusion Eqs. (1) and (2) in the FEM framework can be expressed as a system of linear algebraic equations:

where the matrices ${K}_{x,m},\phantom{\rule{.2em}{0ex}}{C}_{x,m}\left({\mu}_{a}+\frac{i\omega}{c\left(r\right)}\right)$, and F_{x,m} have entries given by:

and the source vector Q_{0} has terms

The source term is defined as a Gaussian distribution, matching the intensity profile at the tip of the optical fiber. Because the source is assumed spherically isotropic, modeling is more accurate when it is centered one scattering distance within the outer boundary. The source vector Q_{m} for fluorescence re-emission is expressed as

and is distributed throughout the domain.

#### 2.2 The inverse model

### 2.2.1 No spatial priors

In the inverse problem, the goal is the recovery of optical properties at each FEM node using surface measurements of light fluence at both the excitation and emission wavelength sequentially (assuming the use of two externally applied sources, one at each wavelength), followed by fluorescence yield reconstruction at the emission wavelength. The computational approach is to minimize the difference between measured fluence, Φ^{Meas}
_{x,m}, at the tissue surface and calculated data, Φ^{C}
_{x,m}, from the model Eqs. (1) and (2) by adjusting the spatial distribution of the unknown parameters through minimization of the ‘objective’ function. The objective function for recovering the optical properties at the excitation wavelength, *μ _{x}* = (

*μ*,

_{ax}*κ*), is given as

_{x}where *NM* is the total number of measurements given by the imaging system, *NN* is the number of parameters representing the optical property distribution which corresponds to the number of nodes in the reconstruction mesh, and *I* is an *NN* × *NN* identity matrix. In general, *χ*
^{2} will not equal zero, but the values of *μ _{x}* for which $\frac{\partial {\chi}^{2}}{\partial {\mu}_{x}}$ is close to zero can be determined, based on an initial estimate of

*μ*. Following the Taylor series method for deriving Newton’s method, $\frac{\partial {\chi}^{2}}{\partial {\mu}_{x}}$ is evaluated at

_{x}*μ*based on an expansion around some nearby point

_{x}*μ*, where the second and higher order terms are ignored, leading to the iterative update equation:

_{x0}where J is the Jacobian matrix, here calculated using the Adjoint method [24]. Equation (13) is known as the Moore-Penrose generalized inverse and is found to be suitable for problems where the number of unknowns to be recovered is much larger than the amount of information (measurements) available. In standard practice, *I* is an identity matrix, and in this work λ is some fixed fraction multiplied by the maximum value on the diagonal of the Hessian matrix *J ^{T} J*, and is therefore updated at each iteration. To recover the optical properties at the emission wavelength,

*μ*= (

_{m}*μ*,

_{am}*κ*), the externally applied illuminating source is changed to one at the emission wavelength and the formulation presented in Eq. (13) is used.

_{m}The recovered optical properties are used in the fluorescence yield reconstruction which also adheres to the minimization formulation presented in Eq. (12). The unknown parameter in Eq. (13) becomes *ημ _{af}* (

*r*), and the Jacobian can be calculated by similar Adjoint properties described above and in Ref. [5].

### 2.2.2 Spatial priors

Spatial prior information is incorporated by assuming a ‘generalized Tikhonov’ penalty term which is similar in structure to Eq. (12) except that the identity matrix is replaced with a Laplacian-type matrix, presented here for the excitation field:

where *NN* is the number of unknowns in the model [16, 25, 26]. The constant, *β* , balances the effect of the parameters with the model-data mismatch in the same manner as λ in Eq. (13). The dimensionless ‘filter’ matrix, *L*, is generated using MRI-derived priors and its construction is flexible. In this application, each node in the FEM mesh is labeled according to the region, or tissue type, with which it is associated (in the MR image). The L-matrix represents a Laplacian-type structure, the diagonal of which is *L _{i,i}*=1 where

*i*is the nodal index. When nodes

*i*and

*j*are in the same region containing

*n*nodes,

*L*=-1/n, otherwise

_{i,j}*L*=0. This effectively relaxes the smoothness constraints at the interface between different tissues, in directions normal to their common boundary. The effect on image quality is similar to that achieved through total variation minimization schemes [27] but easily encodes internal boundary information from MR images. Following the same procedure as in the no-priors case, the parameter update is given by

_{i,j}where much like λ in Eq. (12), *β* is a fixed fraction multiplied by the maximum value on the diagonal of *J ^{T} J*. This update formulation is also used to recover the optical properties at the emission wavelength as well as the fluorescence yield, as described in Section 2.2.1.

### 2.2.3 Reconstruction basis

A number of different strategies for defining the reconstruction basis are possible, including second mesh and pixel basis [28, 29]. The choice of reconstruction basis allows for computational efficiency which serves to reduce the number of unknowns in the algorithm. The problem at hand is twofold: The forward problem requires that the volume of interest is subdivided into adequate number of sub-domains which allow for an accurate description of the calculated fields, whereas a reduction in the number of unknowns improves the ill-posedness of the problem in the reconstruction algorithm. This is addressed by defining a separate reconstruction basis (different from the meshes used in the FEM implementation), upon which the unknowns are updated in Eqs. (13) and (15) [23]. In cases where no prior structural information is available, a pixel basis which defines a set of regularly spaced pixels for the update of the quantities of interest is used. However, when spatial priors are involved, it is crucial to ensure that enough pixels are used to adequately define small structural regions. In this case, a set of regular pixel bases are introduced in each region of interest, which account for the region’s individual shape and size. A semi-adaptive method which allows the number of pixels in each region to be selected was used in the studies presented here. The same region-based pixel basis is used throughout the iterative reconstruction process.

## 3. Methods

#### 3.1 Simulation studies

Test domains for simulation studies were derived from a T1 weighted MR image of a human breast which measured approximately 10.5 cm in diameter. A 2-D slice of the breast volume was discretized into a mesh of approximately 2000 nodes and MR image intensity thresholds were used to assign adipose and fibro-glandular tissue volumes into mesh regions. Figure 1 shows the original MR breast image and its associated discretized, region-labeled mesh. The MR image originated from a clinical exam with our MR coupled frequency domain NIR system and therefore reveals an irregularly shaped breast boundary caused by the fiber optic array slightly compressing the breast at the contact positions. A cancerous tumor region was added as a target anomaly or region of interest for these studies and is depicted in the figure. A second test case with a tumor region located near the center of the domain was used to explore the nonlinear sensitivity of diffuse optical tomographic techniques. The source-detector positions were determined directly from the MR image and represent 16 optical fibers circumscribing the breast domain in a single plane. Light is detected at all non-source fiber positions providing a total of 240 measurements of intensity and phase per wavelength. This simulated configuration matches our experimental fluorescence tomography system.

In this set of studies, regions were assigned values in terms of biologically relevant optical absorbing (HbO, dHb, Water, Lutetium Texaphyrin) as well as scattering parameters (Mie scattering amplitude and power). These values, provided in Table 1, represent typical known in vivo levels of endogenous chromophores and are consistent with previous clinical work [30, 31]. Though it is unknown exactly how Lutetium Texaphyrin (LuTex) will distribute in a human breast, the concentrations used are similar to those reported in ex vivo studies [17] and at least represent a reasonably complex distribution for demonstrating recovery of fluorescence yield. Tissue chromophore concentrations were used to calculate total optical absorption and scattering values at the excitation and emission wavelengths based on experimentally determined values of molar extinction coefficients. Simulated noisy data was generated based on these optical properties in the following manner: 1) frequency domain data for a light source at the excitation wavelength, 2) frequency domain data for a light source at the emission wavelength, and, 3) Continuous Wave (CW) fluorescence emission data at the emission wavelength. This represents a total of 3 data sets for a given imaging session: two frequency domain data sets for determining background optical properties and one CW fluorescence emission data set. CW fluorescence emission data was used to match the capabilities of our experimental system and results in a simplification of Eq. (2) which can be handled by setting *ω* = 0.

Given the modest Stoke’s shift of LuTex, depicted in Fig. 2, the choice of excitation wavelength has practical implications for experimental work. The difficulty in filtering the excitation signal precludes the use of Lutex’s NIR absorption peak (about 735 nm) to excite the fluorophore. In this study, a 690 nm excitation wavelength was used for both simulated and experimental data acquisition. In addition to normally distributed random noise added to the frequency domain data (5% amplitude and 1 degree phase) and CW fluorescence emission (10% intensity), excitation signal leakage through the filter was added to the CW fluorescence emission intensity to simulate a typical 7 OD rejection of the excitation intensity. This number comes directly from experimentally measured rejection estimates for the filters in the tomography system used in the experiments performed here.

The general image reconstruction protocol was as follows:

- Reconstruct for optical properties at the excitation wavelength, μ
_{ax}and μ_{sx}’, with frequency domain data, - Reconstruct for optical properties at the emission wavelength, μ
_{am}and μ_{sm}’, with frequency domain data collected using a laser source at the emission wavelength, - Use the reconstructed optical properties and fluorescence intensity data to recover fluorescence yield.

The same reconstruction algorithm was used to determine background optical properties in steps (1) and (2) and is based on previously reported work [32, 33]. Initial estimates for all parameters were generated using homogenous fitting algorithms which enforce a single value for all nodes. The Jacobian matrix was calculated on a fine mesh of approximately 2000 nodes and interpolated onto a course reconstruction pixel basis for inversion. A 30 by 30 pixel reconstruction basis was used for the no-priors case. The spatial priors reconstruction used a newly developed semi-adaptive pixel basis that redistributes the pixels based on the region information, as described in section 2.2.3. This method ensures that each region contains an adequate number of nodes to approximate the internal structure of the domain. Convergence was defined as less than a 2% change in projection error between successive iterations for the frequency domain optical properties algorithm and less than a 1% change in projection error between iterations for the fluorescence yield reconstruction. Similar algorithmic parameters were used for the phantom studies and are described in further detail below.

#### 3.2 Phantom studies

A spectrometer based tomographic imaging system which couples directly into a Philips 3T MRI magnet was used to acquire fluorescence emission spectra. The system, depicted in Fig. 3, is composed of 16 spectrometers with low noise CCD cameras cooled to -70 C. Sixteen 13 meter long silica fiber bundles circumscribe the imaging domain in contact mode and couple into the spectrometers via custom designed input optics with a six position automated filter wheel. Fluorescence emission signals are processed with long pass interference filters before entering the spectrometer. Each fiber bundle contains eight 400 μm fibers, seven of which directly couple the tissue surface to the spectrometer input while the eighth couples the tissue surface to the light source. This detector configuration minimizes coupling losses and provides parallel detection of full spectra for each source position. Inter-fiber tissue coupling calibration factors were determined prior to imaging using a cylindrical homogenous phantom with a centrally located source/detector. The custom manufactured input optics were designed to optimize filtering and maximize the light detection efficiency of the spectrometers. Camera exposure times can be adjusted to the detected light level to maximize the signal to noise ratio. With 15 detectors per 16 source positions, a total of 240 measurements are recorded for a given acquisition.

Lutetium Texaphyrin provided by Pharamcyclics was diluted in water and used as the imaging fluorophore. The test domain was a solid 5.5 cm diameter hardened epoxy phantom with scatterer and absorbers created by titanium dioxide powder and India ink, respectively [34, 35]. The phantom had a 14 mm hole located approximately 12 mm from the phantom center. The background optical properties were μ_{ax} = 0.005 and μ_{sx}’ = 1.0 mm^{-1}, measured with a frequency domain system near the excitation wavelength. Unlike the simulation experiments, the optical properties were assumed constant throughout the domain in this experiment. These values were also used as the optical properties at the emission wavelength. The hole was filled with a solution of 1% Intralipid to match the scattering value of the background, and varying concentrations of Lutex (0.3125 μM to 5 μM) were added. This represents a simple test case for investigating the imaging response to varying concentrations of fluorophore. The excitation source was a 690 nm laser diode which matches the wavelength used in the simulation studies. Total acquisition time for the fluorescence emission was less than 4 minutes (total of 240 data points).

Even after processing the collected light with a 720 nm long pass interference filter (Omega Optics) which provides 7 OD rejection of the excitation light as well as the filtering offered by the spectrograph grating, emission spectra recorded by the detector are composed of a sum of the pure fluorescence signal and excitation bleed-through. In order to decouple these signals, previously recorded “basis spectra” of the bleed-through signal and the pure fluorescence signal are fit to the data. The process is illustrated in Fig. 4 for one measurement at a single detector. A linear least squares algorithm operating on the basis spectra quantifies the amount of excitation bleed-through and true fluorescence signal that exists in each spectral recording, further reducing the effect of excitation bleed-through. This is accomplished by minimizing the summation

with respect to *a* and *b*, where *y _{i}* is the measured intensity at a given wavelength pixel,

*F*and

*G*are the excitation and fluorescence basis spectra,

*a*and

*b*are the coefficients recovered in the minimization procedure, and

*N*is the number of wavelength pixels per spectrum. Once fit, the fluorescence signal is integrated and becomes the fluorescence emission intensity data for the reconstruction algorithm.

To investigate the improved imaging capability of reconstructing data with spatial priors, each data set was reconstructed with and without spatial “soft” priors. Since the domain was easily characterized geometrically, spatial prior information was determined by direct manual measurement of the phantom. This information was encoded in the fine mesh used for reconstruction. Since background optical properties were previously determined with our frequency domain system, they were assumed known for image reconstruction. Convergence was defined as less than a 1 % change in projection error between iterations. All images were reconstructed using a 2GHz Centrino Duo laptop with 2GB RAM running Windows XP.

## 4. Results

#### 4.1 Simulation results

Figures 5 and 6 display the target values of the two test cases along with images reconstructed using no priors and spatial soft prior information. Qualitatively, image recovery using spatial priors produces significantly more accurate images. Spatial priors preserve the general internal structure of the Lutex distribution, detail that is lost almost entirely in the no-priors case. Images of fluorescence yield show the most dramatic difference between the spatial priors and no priors reconstructions. Without spatial priors, the algorithm appears to have no ability to recover “cancer” regions of elevated fluorescence yield for these complicated cases. However, incorporating spatial priors results in images that qualitatively appear accurate and quantitatively are reasonably close to the true values. Figure 7 provides 1-D cross-sections near the y-axis of each domain. These plots confirm an inability of the no-priors imaging algorithm to recover the simulated fluorescent tumor in either test field. Alternatively, the spatial prior-based imaging algorithm not only picks out the objects of interest, but provides fairly accurate reproductions of the complicated structure of simulated fibro-glandular layers. Mean values for the simulated cancer region near the domain center are 79% and 51% of the true values for the spatial priors and no-priors reconstructions, respectively. These numbers change to 75% and 45% for the case with the cancer region closer to the edge of the domain. They represent a significant improvement in imaging performance; however, they alone do not illustrate the full impact of incorporating spatial priors. The cross-sectional plots indicate that the no-priors images contain virtually no spatial discrimination of the cancer regions. Regional contrasts are depicted in Table 2 and further illustrate a dramatic overall improvement in cancer region quantification with spatial priors.

In addition to improved qualitative and quantitative accuracy, the spatial prior algorithm reduced the reconstruction time significantly. The full reconstruction time for the no-priors case, including background optical property estimation, was just under 9 minutes for both the central and superficial tumor cases. These times were reduced to less than 3 minutes and 90 seconds using spatial priors. In both cases, structural information guided the algorithm to a convergent solution in far fewer iterations than the no-priors cases.

#### 4.2 Phantom results

Figure 8 shows fluorescence yield images recovered from phantom data using both no-prior and spatial prior image reconstruction approaches. A qualitative assessment of the images reveals a dramatic benefit of the spatial prior on image formation. The fluorescent object’s borders are more clearly defined for all concentrations of Lutex when using spatial priors. Furthermore, the values of fluorescence yield throughout the region of interest are more homogeneous, and therefore, more similar to the actual distribution for spatially guided reconstructions. Incorporating spatial priors also suppresses edge artifacts significantly. This is most apparent in images of phantoms with low Lutex concentration. Figure 9 shows a full scale image of the 0.3125 μM phantom. Artifacts near the boundary virtually dominate the no-priors image and represent the largest fluorescent yield values in that imaging domain. None of these artifacts exist in the spatially guided image generated from the same fluorescence emission data. The spatially guided image also includes an easily discernable fluorescent object at the correct location which is difficult to define in the no-priors case. These results indicate that for low fluorophore concentrations in this imaging domain, the correct fluorescent object would be identified only if prior structural information were incorporated in the reconstruction.

Both techniques suffer from a bias in the recovered fluorescence yield which results in a positive value in the region of interest for the case with no fluorophore. This is most likely caused by bleed-through of the excitation light in the emission measurements coupled with an imperfect spectral fitting technique which results in positive values for fluorescence emission intensity even in phantoms containing no Lutex fluorophore. Incorporating additional filtering, such as band-pass source filtering, and optimizing the spectral fitting routine should address the bias signal and decrease the potential for incorrect quantification.

In addition to improving image quality and fluorescence yield quantification, spatial prior-based reconstructions converged in substantially less time than the no-priors algorithms. Reconstruction times ranged from 50 – 90 seconds when not using spatial priors and approximately 22 seconds for spatially guided reconstructions, marking an improvement of just over 75% in some cases. These numbers represent only the fluorescence yield reconstructions and not the recovery of background optical properties since μ_{a} and μ_{s}’ were assumed known from prior measurements. In cases requiring the recovery of background optical properties, the improvement in reconstruction time will be similar to that previously quoted in the simulation results.

## 5. Discussion

This study introduced an effective method to incorporate MR-derived tissue morphology for imaging fluorescence yield at depth in tissue. The conceptual assertion is that diffuse tomography will be more successful as an imaging modality when combined with pre-existing imaging systems which have higher spatial resolution, such as MRI. Simulation and phantom studies were used here to validate the method prior to implementation as a full scale system. Using Lutex phantom data, the soft spatial prior implementation improved qualitative and quantitative imaging response of fluorescence yield. Prior information also suppressed image artifacts and more accurately represented the internal distribution of fluorophore. In the no-priors case, edge artifacts dominate the image for lower concentrations of fluorophore, and provide a misleading interpretation of the internal distribution of the fluorescent agent. The phantom studies represented simple distributions of optical properties and fluorophore concentration and improvements in image quality were still substantial. It is expected that soft prior implementations will benefit imaging performance to an even greater extent in complex tissue domains, an assertion born out in the simulation studies.

A simulated breast mesh derived directly from an MR image of a human breast served as the test bed for a complex tissue domain. The simulations demonstrated no ability to recover the internal distribution of the complicated domain without the use of spatially guided reconstructions. Qualitatively, fluorescence yield images generated without spatial priors had little resemblance to the target domain and completely disregarded the 18mm simulated cancer region, as evidenced by 1-D cross-sectional plots of fluorescence yield. Breakdown of the images in the no-priors case is likely due to poorly recovered background optical properties as well as the complexity of the fluorescence yield distribution itself. Certainly, part of the improvement in fluorescence yield image accuracy can be attributed to improved images of background optical properties.

In these studies, it was assumed that the fluorophore distribution correlates directly to the fatty, fibro-glandular, and tumor tissue layers, which themselves exhibit no significant intra-region heterogeneity, in a given imaging domain. While this assumption appears reasonable based upon the biology of the tissue for endogenous chromophores [30, 31], further studies of the uptake of Lutex in-vivo are required to determine a correlation between tissue type and fluorophore localization. Should this assumption prove to be unrealistic, the “soft” spatial prior approach offers some latitude in terms of correctly identifying structural prior information. The algorithmic implementation groups tissue regions together in a region-specific regularization and allows individual nodes in those regions to update independently. As opposed to “hard” prior approaches where nodal values in a given region are assumed homogenous, a soft prior technique may recover positive fluorescent objects not directly encoded in the spatial prior information. This is a subject of further investigation.

The experimental system introduced here couples directly into a Philips 3T MRI to provide simultaneous MR and NIR fluorescence imaging. Simultaneous imaging simplifies co-registration of the MR image with the NIR reconstruction domains and reduces overall acquisition time. Optimization of this system for in vivo imaging is underway for both small animal and human breast imaging.

## 6. Conclusion

A spatially guided image reconstruction implementation based on prior knowledge of tissue morphology was shown to provide significant improvements in fluorescence yield recovery in complicated tissue volumes and to be highly beneficial for simple domains. Specifically, both phantom and simulation results demonstrated dramatic improvements in recovery and quantification of features in the fluorescence distribution. Structural guidance also reduced image reconstruction time substantially. A newly developed experimental system couples full-volume deep-tissue fluorescence spectroscopy capabilities directly into the MR bore for simultaneous MR-NIR fluorescence image acquisition. The results presented here show promise for this approach in all cases and tissue volumes considered. Extensions of this study are underway to determine the effect of incorrectly identifying the structural prior, especially in cases where the MR images produce false negative or false positive readings. A full characterization of the imaging limits of the experimental system will further complement the results presented here. It is clear that incorporating anatomical features derived from MR images in DOFT image reconstruction will improve sensitivity to lower concentrations of fluorophore, qualitative accuracy, and fluorescence yield quantification in-vivo.

## Acknowledgments

This work was funded by the National Institutes of Health grants RO1 CA109558, RO1 CA69544, U54 CA105480, as well as Department of Defense Breast Cancer pre-doctoral fellowship BC051058.

## References

**1. **D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yeild and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. **36**,2260–2272 (1997). [CrossRef] [PubMed]

**2. **H. B. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite- element-based algorithm and simulations,” Appl. Opt. **37**,5337–5343 (1998). [CrossRef]

**3. **D. J. Hawrysz and E. M. Sevick-Muraca, “Developments toward diagnostic breast cancer imaging using near-infrared optical measurements and fluorescent contrast agents,” Neoplasia **2**,388–417 (2000). [CrossRef]

**4. **M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three-dimensional, Bayesian image reconstruction from sparse and noisy data sets: near-infrared fluorescence tomography,” Proc. Natl. Acad. Sci. USA **99**,9619–9624 (2002). [CrossRef] [PubMed]

**5. **A. B. Milstein, O. Seungseok, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomograhy,” Appl. Opt. **42**,3081–3094 (2003). [CrossRef] [PubMed]

**6. **V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. **26**,893–895 (2001). [CrossRef]

**7. **V. Ntziachristos and R. Weissleder, “Charge-coupled-device based scanner for tomography of fluorescent near-infrared probes in turbid media,” Med. Phys. **29**,803–809 (2002). [CrossRef] [PubMed]

**8. **E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. **30**,901–911 (2003). [CrossRef] [PubMed]

**9. **A. Godavarty, A. B. Thompson, R. Roy, M. Gurfinkel, M. J. Eppstein, C. Zhang, and E. M. Sevick-Muraca, “Diagnostic imaging of breast cancer using fluorescence-enhanced optical tomography: phantom studies,” J. Biomed. Opt. **9**,488–496 (2004). [CrossRef] [PubMed]

**10. **R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissue with noncontact measurements,” IEEE Trans. Med. Imaging **23**,492–500 (2004). [CrossRef] [PubMed]

**11. **S. V. Patwardhan, S. R. Bloch, S. Achilefu, and J. P. Culver, “Time-dependent whole-body fluorescence tomography of probe bio-distribution in mice,” Opt. Exp. **13**,2564–2577 (2005). [CrossRef]

**12. **A. Godavarty, M. J. Eppstein, C. Zhang, and E. M. Sevick-Muraca, “Detection of single and multiple targets in tissue phantoms with fluorescence-enhanced optical imaging: feasibility study,” Radiol. **235**,148–154 (2005). [CrossRef]

**13. **B. Brooksby, S. Jiang, C. Kogel, M. Doyley, H. Dehghani, J. B. Weaver, S. P. Poplack, B. W. Pogue, and K. D. Paulsen, “Magnetic resonance-guided near-infrared tomography of the breast,” Rev. Sci. Instrum. **75**,5262–5270 (2004). [CrossRef]

**14. **X. Intes, C. Maloux, M. Guven, B. Yazici, and B. Chance, “Diffuse Optical tomography with physiological and spatial a priori constraints,” Phys. Med. Biol. **49**,N155–N163 (2004). [CrossRef] [PubMed]

**15. **M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. **50**,2837–2858 (2005). [CrossRef] [PubMed]

**16. **P. K. Yalavarthy, H. Dehghani, B. W. Pogue, C. M. Carpenter, H. B. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” IEEE Trans. Med. Imaging **In review** (2006).

**17. **G. Kostenich, A. Orenstein, L. Roitman, Z. Malik, and B. Ehrenberg, “In vivo photodynamic therapy with the new near-IR absorbing water soluble photosensitizer lutetium texaphyrin and a high intensity pulsed light delivery system,” Photochem. & Photobiol. **39**,36–42 (1997). [CrossRef]

**18. **K. W. Woodburn, Q. Fan, D. R. Miles, D. Kessel, Y. Luo, and S. W. Young, “Localization and efficacy analysis of the phototherapeutic lutetium texaphyrin (PCI-0123) in the murine EMT6 sarcoma model,” Photochem. & Photobiol. **65**,410–415 (1997). [CrossRef] [PubMed]

**19. **M. Zellweger, A. Radu, P. Monnier, H. van den Bergh, and G. Wagnieres, “Fluorescence pharmacokinetics of Lutetium Texaphyrin (PCI-0123, Lu-Tex) in the skin and in healthy and tumoral hamster cheek-pouch mucosa,” Photochem. & Photobiol. B **55**,56–62 (2000). [CrossRef]

**20. **A. Synytsya, V. Kral, P. Matejka, P. Pouckova, K. Volka, and J. L. Sessler, “Biodistribution assessment of a lutetium(III) texaphyrin analogue in tumor-bearing mice using NIR Fourier-transform Raman spectroscopy,” Photochem. & Photobiol. **79**,453–460 (2004). [CrossRef] [PubMed]

**21. **H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and P. K. D., “The effects of internal refractive index variation in near infrared optical tomography: A finite element modeling approach,” Phys. Med. Biol. **48**,2713–2727 (2003). [CrossRef] [PubMed]

**22. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. **20**,299–309 (1993). [CrossRef] [PubMed]

**23. **K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. **22**,691–701 (1995). [CrossRef] [PubMed]

**24. **S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part2: Finite-element-method calculations,” Appl. Opt. **34**,8026–8037 (1995). [CrossRef] [PubMed]

**25. **A. Borsic, W. R. B. Lionheart, and C. N. McLeod, “Generation of anisotropic-smoothness regularization filters for EIT,” IEEE Trans. Med. Imaging **21**,579–587 (2002). [CrossRef] [PubMed]

**26. **B. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. STQE **9**,199–209 (2003).

**27. **K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,” Appl. Opt. **35**,3447–3458 (1996). [CrossRef] [PubMed]

**28. **K. D. Paulsen, P. Meaney, M. Moskowitz, and J. Sullican Jr., “A dual mesh for finite element based reconstruction algorithms,” IEEE Trans. Med. Imaging **14**,504–514 (1995). [CrossRef] [PubMed]

**29. **M. Schweiger and S. R. Arridge, “Optical tomographic reconstruction in a complex head model using a priori boundary information,” Phys. Med. Biol. **44**,2703–2722 (1999). [CrossRef] [PubMed]

**30. **B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. USA **103**,8828–8833 (2006). [CrossRef] [PubMed]

**31. **S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J. Gibson, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “In vivo hemoglobin and water concentrations, oxygen saturation, and scattering estimates from near-infrared breast tomography using spectral reconstruction,” Acad. Radiol. **13**,195–202 (2006). [CrossRef] [PubMed]

**32. **T. O. McBride, B. W. Pogue, S. Jiang, U. L. Osterberg, K. D. Paulsen, and S. P. Poplack, “Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging,” Opt. Lett. **26**,822–824 (2001). [CrossRef]

**33. **B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. **38**,2950–2961 (1999). [CrossRef]

**34. **S. Jiang, B. W. Pogue, T. O. McBride, and K. D. Paulsen, “Quantitative analysis of near-infrared tomography: sensitivity to the tissue-simulating precalibration phantom,” J. Biomed. Opt. **8**,308–315 (2003). [CrossRef] [PubMed]

**35. **B. Pogue and M. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. **11**,0411021–04110216 (2006). [CrossRef]