Near-infrared (NIR) biomedical optical imaging depends upon recovering an interior optical property map from measurements conducted at the air–tissue interface in order to image diseased tissue volumes. In the past decade, several investigators have sought methods to measure NIR propagation and to recover optical property map mathematically to differentiate diseased tissues from normal tissues based upon the optical contrast due to absorption and isotropic scattering. 1 The ability to optically image or to detect diseased tissue volumes located deep within tissues depends upon the endogenous contrast provided by differences in tissue absorption (owing primarily to hemoglobin) and in scattering (presumably owing to changes in cell density or size).
However, using in vitro measurements Troy et al.; 2 showed that the optical contrast due to normal and diseased human breast tissues may not be sufficient for the successful detection of breast cancer. In contrast, recent two-dimensional (2D), quantitative image reconstruction of in vivo absorption within human breast indicate that absorption contrast between normal and malignant tissues may be twofold and indeed sufficient for breast cancer screening. 3 The in vitro results from Fishkin et al.; 4 and Fantini et al.; 5 also suggested that the optical properties of malignant and normal tissues may be dramatically different. Indeed the spectroscopy results and the optical mammograms reported by Franceschini, 6 Moesta, 7 and Grosenick 8 suggest that breast tumors can be detected by endogenous absorption contrast. Yet these positive results are obtained predominantly with palpable legions and/or lesions positively identified by conventional x-ray mammography. Recently, many groups have engaged in efforts to detect smaller lesions and to enhance the optical contrast by the use of exogenous contrast agents. Sevick et al.; 9 and Li et al.; 10 showed that the exogenous contrast offered by fluorescent compounds is superior to that provided by nonfluorescing, light-absorbing compounds when time-dependent photon migration measurements are employed. While the preferential uptake of fluorescent contrast agents into disease tissue volume of interest can be primarily responsible for contrast, the kinetics of fluorescence decay processes can be environmentally specific to tissue volumes and can further induce optical contrast for the detection. To date, the reconstruction of internal fluorescent properties of quantum efficiency and lifetime from synthetic data has been difficult, especially when the finite partitioning of fluorescent compounds takes place between simulated normal and diseased tissues. The distribution of background fluorophore in the simulated “normal” tissues tends to serve as many secondary emission sources, thus contributing to the total detected signal and making it more difficult to differentiate the location and size of a simulated “heterogeneity.” Not surprisingly, there has been no demonstrated reconstruction from synthetic data without a priori information or estimates of correct absorption coefficient map for successful recovery of lifetime. 11 12 13 In addition, while the added contrast is promised from fluorescent agents, there has been little or no attention to the offsetting reduction in signal to noise ratio (SNR) of fluorescence signals in comparison to signals measured at the wavelength of the incident light. More importantly, since one can expect a large dynamic range of fluorescence signals associated with sources and detectors located near and far from a fluorophore laden “heterogeneity,” issues of noise tolerance of the fluorescence-enhanced absorption inversion algorithm represent major challenges that have only recently begun to be addressed. 14
In this paper, we restrict our attention to the recovery of absorption cross section images from fluorescence frequency- domain photon migration (FDPM) signals arising from a fluorophore whose lifetime is assumed constant, such as Indocyanine Green (ICG)—an agent approved by the Food and Drug Administration for use in assessing hepatic function and retinopathy. We term this imaging, “fluorescence-enhanced absorption imaging.” We first assess the experimental errors of single-pixel FDPM measurements in a homogenous phantom at the incident and emission wavelengths as a function of signal level. Next, we represent the experimental errors as a SNR, generate synthetic measurements corrupted by the representative experiment error, and use these data sets to test the noise tolerance of our inversion algorithms. In our inversions, the unknown parameters governing the generation and propagation of fluorescent light are the absorption and scattering coefficients at excitation and emission wavelengths μa x,μa m and μs x′ ,μs m′ and the product of the absorption coefficient owing to fluorophore and its quantum efficiency ϕμax→m.
Our contribution is organized as follows: First, the formulation of the reconstruction algorithm for the recovery of absorption from FDPM measurements at the incident wavelength and the emission wavelength. The experimental assessment of FDPM measurement error at excitation and emission wavelengths is described in Sec. 2. In Sec. 3, the measurement setup for the experimental assessment of SNR of FDPM and the procedures for the forward and inverse problem with synthetic data sets are then described. In Sec. 4, we present our assessment of SNR of FDPM and then present the corresponding reconstruction results as a function of SNR for synthetic FDPM data sets at the incident and emission wavelengths. Finally, in Sec. 5 we conclude with a “laundry list” for the challenges and potential solutions for fluorescence-enhanced absorption imaging.
In this section, the formulation of inverse algorithms to recover the optical property μa x at incident and emission wavelength is briefly discussed. Since the goal of this investigation is to assess the measurement noise tolerance of inverse algorithms, a simple framework for the experimental assessment of the SNR of the detection system is presented.
Reconstruction of μ a x from FDPM Measurements at Incident Wavelength
In highly scattering media such as tissues, the propagation of light is well described by the diffusion approximation to the radiative transport equation. 15 The excitation photon fluence Фx(r,ω) at angular frequency ω, and detected at position r is described by the diffusion equation in the frequency domain for excitation light (superscript x):μa x(r) is the absorption coefficient, Dx(r)=[3⋅(μa x(r)+μs x′ (r))]−1 is the optical diffusion coefficient, c is the speed of light in the medium, Sx(r,ω) is the excitation source term, and μs ′x(r) is the isotropic scattering coefficient. kx 2(μa,μs)=3(μa+μs)[−μa+(iω/c)] and S=3(μa+μs)δ(r−r s). Here, the gradient of the optical diffusion coefficient is considered negligible—a realistic assumption for biomedical imaging given that the isotropic scattering does not alter dramatically between normal and diseased tissues.
Reconstruction of μ ax → m from FDPM Measurements at Emission Wavelength
The complex emission light fluence Фm(r,ω) is described similarly to the excitation fluence with the exception that the source of emission photons is distributed within the volume and coupled to the excitation fluence Фx(r,ω):μax→m is the absorption coefficient due to the fluorophore; Sm(r,ω) is the emission fluence source term; Dm(r)=1/[3(μa m+μs m′ )] is the diffusion coefficient at the emission wavelength λm, and τ(r) and is the probe lifetime at position r. Since we assume the lifetime of ICG does not change from the background value upon uptake to the diseased tissue, we will denote the constant lifetime as τ0 henceforth. Also, we assume the absorption coefficient due to fluorophore μax→m is directly proportional to the concentration of fluorophore Nt. For the purpose of our investigation, we assume that the dye’s fluorescence and absorption spectra are well separated so that we can ignore the small influence of the excitation of fluorophores by the fluorescent remission. In addition, first order single exponential fluorescence decay kinetics is assumed.
Since the emission diffusion equation is in the form of an inhomogeneous differential equation, the Green’s function is used to obtain the analytical solution of the emission fluence Фm(r). Equation (3) can be expressed as follows:[∇Dm(r)/Dm(r)]⋅∇Фm(r) in Eq. (4) accounts for the discontinuity in Dm(r). However, since there is little or no variation of Dm(r) in the emission wavelength λm, this term is negligible. The Green’s function corresponding to Eq. (4) consequently satisfies 4) and (6), and with the use of Green’s theorem, the emission fluence Фm(r d,r s) at the detector position r d, resulting from an excitation source at position r s, can be expressed as follows: 7) one obtains Ф˜m=FX, can be denoted as F∈CM×N, X∈R N, Ф˜m∈CM, respectively, and M=KL.
Experimental Determination of the SNR
To evaluate the performance of the inverse algorithm, the image reconstruction is usually carried out using a synthetic data set generated by the forward simulator. It is common practice to add random noise to the synthetic forward solution to mimic the experimental noise. The framework for the measurement of SNR of FDPM is presented in this section. Experimentally assessed SNR at both excitation and emission wavelengths are later used as the guideline for the level of random noise added to the synthetic data sets for the reconstruction algorithms.
The actual signal that is detected at the air–tissue interface and used as input to the inversion algorithm is proportional to the gradient of the complex fluence at either incident or emission wavelengths ∇Ф(r d,ω) plus a contribution owing to measurement error. When zero-fluence or partial current boundary conditions are employed, the gradient of the complex fluence is also directly proportional to the complex fluence itself.
Experimental values of ac amplitude and phase delay can be represented as a complex signal S=M AC exp(iθ)=X+iY. With a large number of FDPM measurements, the sample means of the real part and imaginary part of S were used as the true signal levelX¯ and Y¯ are the sample mean of real and imaginary parts, respectively, and Nx and Ny represent the noise in the real and imaginary parts. The expectation value of the signal is then simply: σX 2 is the variance of X and σY 2 is that of Y. We can then calculate SNR to be
To determine the SNR from experimental FDPM measurements, the variances and means of phase and ac amplitude were determined from repeated experimental measurements and then used to compute the signal power as well as the variances of the real and imaginary parts of the fluence σx 2 and σy 2.
Methods and Approaches
In this section, the experimental setup to measure SNR at incident and emission wavelengths is briefly described. Also, the simulated phantom and the generation of synthetic data and the inversion strategy are outlined.
Experimental Setup for SNR Measurements at Incident and Emission Wavelengths
To validate the noise model we used in the previous sections, frequency domain measurements are conducted in a homogeneous scattering medium which both fluoresced and absorbed. Figure 1 illustrates the schematic diagram of the measurement equipment. A laser diode (70 mW at 785 nm, Sanyo) was intensity modulated at a rf of 100 MHz with a frequency synthesizer (signal generator 2002D, Marconi Instruments). The beam splitter was used to direct 10 of the incident beam to the reference photomultiplier tube (PMT) (Hamamatsu R928). Both sample and reference beams were collimated and delivered to 1000 μm core optical fibers (Thor Lab 0.39-NA TECS multimode fiber). PMTs were gain modulated at the modulation frequency of a laser diode plus an offset frequency of 100 Hz. The heterodyned PMT signals were analyzed for phase shift and ac and dc amplitude using PC-based LABVIEW software (National Instruments Corp., Austin, TX). A continuously variable neutral density filter wheel (Newport, Irvine, CA) was used to change the dc intensity of the source. The source and detector fibers were positioned at the side of an 8×8×8 (cm 3 ) phantom separated by 1 cm, which was filled with 1 Intralipid solution. For FDPM fluorescence measurements, an Intralipid solution containing ICG at a concentration of 0.1 μM was employed. The absorption coefficient and the isotropic scattering coefficient are 0.0029 and 10.75 cm−1, respectively. A narrow band interference filter at 830 nm (10 nm full width at half maximum) was positioned inside the sample PMT housing to collect the light at the emission wavelength.
The detected FDPM signals as a function of incident dc level were sampled 100 times with a 0.2 s scanning time. dc, ac amplitude, and phase shift relative to the source were measured and their means and standard deviations computed. If the detector noise is limited by shot noise statistics, the variances of dc and ac will increase with dc amplitude. 18 We also compute the signal to noise ratio of both excitation and emission measurements.
Simulated Phantom for Forward and Inverse Problems
The dimension of the 2D square phantom for the forward and inverse problems was chosen to be 4 cm×4 cm discretized into a 33×33 grid. This leads to the distance between neighboring pixels of 0.125 cm. Since the source located at the air–tissue interface can be modeled as a single isotropic source at one transport mean free path from the boundary, 19 we modeled the source position one pixel length from the medium boundary. Fourteen sources and detectors are represented around periphery of the phantom at equivalent separation distances. The refractive index of the phantom was set to 1.4 to mimic the properties of tissue. To account for the refractive index mismatch at the air–tissue interface, a partial-current boundary condition is used instead of a mathematically simpler zero-boundary condition. 19 20 In this work, we consider a realistic 10:1 uptake ratio of fluorescent agent leading to ten times greater absorption within the heterogeneity as in its surrounding “background.” We report results for 100 MHz excitation and assume no fluorescent lifetime changes between the heterogeneity and its surrounding “background.” The heterogeneity is 1 cm2 and is positioned at (18, 18) in the 33×33 grid. We assume constant scattering properties throughout the phantom since independent work by others show unmodeled variations in scattering have little influence on reconstructions. 21
Generation of Synthetic Data Sets and Solution of the Inverse Problem
For the solution of the forward problem, the known optical properties of background and the heterogeneity shown in Table 1 are used to calculate the fluence at each node, and the complex fluence values are collected at the detector nodes. Gaussian random noise is added and the calculated values are used as simulated measurement input into the inverse algorithm. IMSL FORTRAN routines were used to add the random noise to the simulated measurement
|Optical properties of background and heterogeneity in the simulation.|
N(0,1) represents a random number of Gaussian distribution between 0 and 1.
For the reconstruction of μax(r) at the incident wavelength, the BIM is used. In order to reconstruct the spatial map of μax→m(r) detailing the heterogeneity, we discretize Eq. (7) into a series of equations in terms of Gf, Фx, and measurements of Ф˜m(r d,r s). The solution of the forward problem and determination of Gf and Фx is obtained using MUDPACK, 22 a finite difference differential equation solver for FORTRAN. Both forward and inverse problems are carried out in a two-dimensional geometry, since the three-dimensional geometry is prohibitive for the preliminary computations presented herein. For the simulated measurements, we consider the excitation source to be amplitude modulated by a frequency of ω. Measurements of phase shift θm and amplitude of ac component Mm [or Фm(r d,r s)=Mmeiθm ] are obtained at detector position r d in response to excitation source at r s.
Using the synthetic data sets of FDPM measurements at both the incident and emission wavelengths, the spatial map of the absorption coefficient due to fluorophore μax→m(r) is obtained by the Levenberg–Marquadt method. 6 23 The update Δμax→m is obtained by solvingF(i)F(i)H matrix since the matrix is ill posed. 24 The superscript H denotes a complex conjugate. To monitor the convergence of the inversion algorithm, the root-mean-square-error (RMSE) of phase, defined as θ meas (j) is the measured phase of the fluence from the jth source/detector pair and θc,n(j) is the computed fluence at nth iteration due to the jth source/detector pair. RMSE of amplitude RMSE M μa x or μax→m
Results and Discussion
Experimental Assessment of FDPM Measurement Error at Incident and Emission Wavelengths
Figures 2 and 3 illustrate the computed values of SNR versus the dc power of the source. For FDPM measurements at the incident wavelength, the SNR was computed to be relatively constant around 50 dB. For FDPM measurements at the emission wavelength, the SNR level was around 35 dB. The reduction in SNR can be attributed to the weaker signal level of emission light and the increased gain setting of the sample PMT required to conduct the measurement. Nonetheless, the SNR deteriorates with decreasing signal power, indicating that it might be important to consider variable measurement error within a robust strategy for image recovery in fluorescence-enhanced imaging. For the purposes of this work, we assume a constant SNR for FDPM measurements at both incident and emission wavelengths to test the noise tolerance of fluorescence-enhanced absorption imaging. Nonetheless, this data shows that Bayesian reconstruction approaches, which account for variable measurement error such as those developed by Eppstein et al.;, 14 may be appropriate for fluorescence-enhanced absorption imaging.
Reconstruction of μ ax from the Scattered Incident Wave
Using the measured SNR in the incident and emission wavelength, optical properties were reconstructed to test the noise tolerance of the inversion algorithms. The procedures to add experimentally relevant noise and the characterization of noise in the experiment are discussed in the Appendix.
Figure 4 shows the reconstruction results of absorption coefficient from FDPM synthetic data at the incident wavelength. The initial guess of homogeneous background optical properties was employed. The reconstruction algorithm is able to recover the optical map with SNR of 40 dB [Fig. 4(b)] and 50 dB [Fig. 4(c)], but as the SNR becomes smaller, the reconstruction produces unsatisfactory results when the SNR level is below 35 dB [Fig. 4(a)]. At a SNR level of 35 dB and below, the reconstruction failed. Even with a perfect absorber in the lossless media, the dominant first-order scattered wave, Ф scatt x(r d), is small in comparison to the incident wave Ф inc x(r d). Hence, the contrast due to the absorbing heterogeneity is expected to be small. 9 As the optical property between the heterogeneity and its surrounding decreases, the magnitude of the scattered wave becomes even smaller. When the heterogeneity is buried deep within the tissue, the contrast available for the successful image reconstruction diminishes even more and can be lost in the noise of phase and amplitude measurements. Furthermore, due to the ill-posedness of the problem, the reconstruction results becomes spurious when SNR is smaller than 35 dB in spite of regularization. When the regularization parameter is increased, the pseudoinverse of the Jacobian matrix becomes too diagonally dominant useless results are produced. On the other hand, if the regularization parameter is set too small, the additive noise leads to spurious reconstruction.
Figure 5 illustrates plots of iteration versus RMSE values associated with the absorption maps. Both RMSE θ and RMSE M illustrate convergence after ten iterations while RMSE μa shows that increased accuracy, of course, associated with the highest fidelity data sets. When SNR levels fall below 35 dB, both RMSE θ and RMSE M values increase after 2–3 iterations, leading to unsuccessful reconstruction of the absorption coefficient distribution.
Reconstruction of μ ax → m from the Emission Wave
Figure 6 shows that the fluorescence-enhanced absorption imaging algorithm was able to recover the absorption map owing to the absorption contrast from fluorophores, μax→m, when SNR varied from 50 to 20 dB. Figure 6(a) shows that even with the low SNR of 20 dB, the reconstruction of μax→m was possible. This robustness of algorithm is attributable to the added phase and amplitude contrast from the fluorescence decay kinetics, which is absent from the conventional absorption algorithm. The greater the partitioning of the fluorescent dye within the heterogeneity, the stronger will be its fluorescent source enabling optimal contrast. Fluorescent dyes, whose quantum yield and lifetime changes upon the preferential uptake by the diseased tissue, can further aid contrast. Even if the quantum efficiency ϕ and lifetime τ do not change in the heterogeneity and in its surroundings, the ratio of μax→m and the additional phase shift due to lifetime τ will give rise to the larger contrast for measurement at the fluorescent wavelengths than at the incident wavelength.
The influence of the added noise is seen in Fig. 7, which illustrates the RMSE θ versus iteration plots [Fig. 7(a)]. Unlike the absorption imaging cases where the reconstruction failed below SNR of 35 dB, the fluorescence-enhanced imaging algorithm recovered the optical map even at the SNR of 20 dB. Our findings concur with Sevick-Muraca et al.; 9 and Li et al.; 10 in their investigation of the detection limit of fluorescent inhomogeneities in turbid media. They concluded that, for a given fluorophore concentration and object size, both amplitude and phase perturbation of emission wave are greater than those measured at the incident wavelength.
Summary and Conclusion
While the fluorescence-enhanced absorption imaging algorithm presented herein is similar in structure to the Born iterative method, it nonetheless differs in that it employs FDPM measurements at the emission rather than at the incident wavelengths. In the past, Sevick et al.; 9 and Li et al.; 10 have shown that the additional phase and amplitude information associated with the generation of a fluorescent wave from NIR excitable exogenous contrast agents can enhance contrast and can aid in the reconstruction of optical property maps of the diseased tissues. In this work, we furthermore show that fluorescence-enhanced contrast is advantageous even when the tradeoff of reduced SNR is considered.
Using a typical single-pixel FDPM system to interrogate a homogeneous phantom at differing levels of source power, we found that the reduction of SNR in emission FDPM measurements was not great enough to discount the extra contrast owing to fluorescence decay kinetics. Furthermore, we found that while SNR was relatively constant for FDPM measurements made at the incident wavelength, it varied with signal power level at emission wavelengths. While SNR levels can be expected to vary depending upon FDPM instrument design (such as the case in multipixel FDPM instruments), 25 they nonetheless should be measured and addressed within the framework of the inverse imaging algorithm used to map the pertinent tissue optical properties. Furthermore, in cases where there is a wide dynamic range of detected signal levels (as in the case of detection of the incident wave at a detector close to a large absorbing heterogeneity or in the case of the detection of an emission wave far from a fluorescently enhanced tissue heterogeneity) SNR may not be expected to be constant. In these cases, the incorporation of noisy FDPM measurements into an inverse algorithm without appropriate weighting against low SNR will limit the number of measurements and consequently the amount of information for image recovery. Our results point to the continued development in biomedical optical imaging of: (i) NIR excitable contrast agents; (ii) high fidelity data collection; (iii) quantification and assessment of SNR; and (iv) Bayesian inversion strategies which provide dynamic regularization based upon assessed measurement error.
Appendix: Characterization of Noise in the Experiment
Suppose Ф signal (r) is the “true” fluence at the position r and N(r) is the independent and additive random noise normally distributed. If it is assumed that the measurement is shot-noise limited, then independent and random Gaussian noise is added to the complex fluence which is sufficient to represent measurements in a synthetic data set. Consequently, the synthetic fluence is the sum of Ф signal (r) and N(r)Ф signal (r) is the complex number quantity represented by ac amplitude M and by the phase delay θ, N(r) should also have real and imaginary components N(r) |M|2=[abs(Ф signal )]2, σnr 2, and σni 2 are variances of the real and imaginary part of N(r), respectively. We use Eq. (A4) to generate synthetic data sets as well as to determine the SNR from variances of phase and ac measurements.
Figures 8(a) and 8(b) show the variance of dc and ac as a function of detected dc when FDPM measurements were conducted at the incident wavelengths. Figures 8(c) and 8(d) provide the corresponding measurements conducted at the emission wavelength. In both cases, variances in dc and ac amplitude increased as the dc power of the source increased. These plots indicate that the detection system is shot noise limited. The noise in the real and imaginary parts are relatively independent with a correlation coefficient of ∼0.1 for the measurements made at the incident wavelength and ∼0.05 for those at the emission wavelength. These results justify our assumptions made in the generation of synthetic data sets.