1. Introduction

Fluorescence-enhanced optical imaging may hold promise for diagnostic cancer imaging with molecularly targeting fluorescent agents. Currently, the “gold-standards” for clinical molecular imaging are the nuclear techniques of gamma scintigraphy, positron emission tomography (PET), and single photon emission computed tomography (SPECT) [1, 2]. Since radiotracers have a finite half-life and cannot be “re-activated” in vivo, the signal to noise available for planar imaging or tomographic reconstruction can be limiting. A similar scenario exists for bioluminescence imaging which can also be starved for signal since the generation of light is governed by the diffusional encounter of an enzyme and its consumable substrate. In contrast, fluorescence imaging may (i) have higher signal to noise ratio owing to the ability of a fluorophore to undergo repeated activation and radiative relaxation; (ii) result in enhanced target-to-background ratio (TBR) values since imaging can be delayed without concern of agent half-life or availability of substrate; but (iii) may suffer greater penetration losses at increasing tissue depths than nuclear imaging approaches [3]. The opportunity to conjugate fluorophores to targeting peptides, antibodies, and enzyme substrates has already been well demonstrated in several laboratories using small animal imaging.

When fluorophore excitation and emission occurs in the near-infrared range (750-900 nm) and light is multiply scattered as it travels a centimeter or more in tissues, image formation can be achieved from inversion of time-dependent measurements made at the tissue surface using the coupled diffusion equations [3–7]. When used in small volumes, tomographic imaging requires the use of the time-dependent radiative transport equations (RTE) [8], although some investigators have employed diffusion based tomography with time-invariant measurements [9]. In addition, while diffusion-based tomography is performed across large volumes with measurements between points of illumination and collection on the tissue boundary, 3-D image reconstruction has also been developed for area illumination and collection [10, 11]. Herein, we will focus on diffusion-based tomographic imaging from time-dependent measurements conducted in the frequency-domain between points of illumination and collection in the presence of natural anatomical heterogeneity and heterogeneous distribution of fluorescent agent in the geometry typical of human breast. We chose this geometry owing to our past experimental studies involving a breast-shaped phantom [12] and employ optical property values that are similar to those reported in the literature [13–23] and summarized in Table 1.

In the following, we present in the Methods section: (i) the forward model and finite element mesh for predicting boundary measurements of frequency-domain photon migration (FDPM) measurements between points of illumination and collection in the presence of optical property heterogeneity, and (ii) a brief description of the inversion scheme used to provide reconstructed images. In the Results and Discussion section, we provide example reconstructed images at varying levels of optical property heterogeneity to show the insensitivity to endogenous optical properties in fluorescence-enhanced optical imaging. Our results also indicate the sensitivity of tomographic imaging to uneven distribution of molecularly targeted agents in normal tissues. The results provide the framework for OAIQ.

2. Methods

Time-dependent, frequency-domain measurements consist of launching intensity modulated excitation light (typically modulated at 100 MHz) from a single point on the tissue boundary. The intensity modulated excitation light propagates as a photon density wave through the tissue, activating fluorophores and generating intensity modulated emission light. The generated emission photon density wave is phase-shifted and amplitude attenuated relative to its activating excitation light due to the decay kinetics of the fluorophore. The emission photon density wave propagates to the tissue boundary where it is collected for evaluation of its amplitude and phase delay for input into the image reconstruction algorithm.

2.1 Forward model and finite element solver

Near-infrared light propagation in tissues can be modeled by the diffusion approximation of the radiative transport equation. In the frequency domain the photon diffusion equation is generally written as:

Table 1. Experimental breast optical property values reported in literature

View Table | View all tables in this article

Here ω is the modulation frequency of the NIR source (rad/s); Φ(r, ω) is the photon fluence rate (photons/(cm²s)) at position r ; D is the photon diffusion coefficient (cm); μ_a is the absorption coefficient (cm ^-1); S(r,ω)is the photon source strength (photons/(cm ³ s)) at position r; c is the speed of light in the medium (cm/s). The generation and propagation of the fluorescence diffuse photon density wave can be described by the following coupled diffusion equations:

Where the subscript x denotes excitation and m denotes emission. The term μ_ax,mi denotes the absorption due to endogenous chromophores; μ_ax,mf denotes the absorption due to the exogenous fluorophores; ϕ is fluorescent quantum; and τ is fluorescent lifetime (s). The diffusion coefficient is given by:

Where μ_sx,m denotes the scattering coefficients at excitation and emission wavelength; and g is the coefficient of anisotropy of the medium. We solve these equations with the Robin type boundary conditions:

Where n denotes the outward normal to the surface and γ is the constant depending upon the optical refractive index mismatch at the boundary. Eq. (2), (3), (4) and (5) can be solved numerically to yield,

Here Φ_x,m is a complex number, and θ_x,m is the measured phase lag and I_ACx,m is the measured amplitude of the photon density wave at excitation and emission wavelengths.

Unlike finite difference methods, finite element methods can solve coupled diffusion equations Eq. (2) and (3) over complex domains. We employ a Galerkin scheme for solution of these equations over a breast-shaped geometry described in Fig. 1 and used in prior experimental studies in our laboratory [24].

2.2 Endogenous and exogenous optical property heterogeneity

We consider the lumpy-object model developed by Rolland and Barrett [25] to simulate the normal background anatomy as a representation of the non-specific distribution of the fluorescent agent as well as the natural heterogeneity of the endogenous tissue optical properties. The lumpy backgrounds consist of a random number of single structures or “blobs” located at random locations. Mathematically these structures can be represented as:

Where b(r) is lumpy background; b ₀ is the spatial mean of the background; N_p is the Poisson-distributed number of lumps; and r _n is the uniformly distributed location of nth lump. The function lump has the following form:

Where l₀ is the lump strength, w is the lump width, Ω is the domain, and V(Ω) is the volume of the domain. The second term in this expression satisfies the requirement that the mean of the lumpy background is equivalent to b ₀, i.e., 〈b(r)〉 = b ₀ [26].

Lumpy backgrounds have the advantage of being mathematically/statistically tractable and stationary [27]. This facilitates an organized manner for implementing OAIQ tools for image analysis. These lumpy backgrounds are represented in the uniform mesh of the FEM forward model of generated image data sets which are used to compute the tomographic reconstruction images.

 

Fig. 1. The geometry of the breast-shaped phantom. The dimensions are in centimeters. The bigger cylindrical volume has a circular base of diameter 20 cm. The hemisphere has a radius of 5 cm.

Download Full Size | PDF

The breast geometry (shown in Fig. 1) consisted of two cylinders of radii 10 cm and 5 cm with heights of 2.5 cm and 0.5 cm respectively, and a hemisphere of radius 5 cm. A total of 6956 nodes and 34413 tetrahedral elements were used to discretize the domain. The lump strength, l₀, was set to be a prescribed percentage ranging from 0 to 100% of the mean background value b₀; and the lump spread, w, and number of lumps, N_p, were set to be 5 mm and 100 respectively. As an example, 5% lumpy background of the endogenous or exogenous properties meant that the strength of lump, l₀, was 0.05 times the average background value (b₀) of the optical property. Due to the large volume and surface illumination of excitation light on the hemispherical portion, the intensity of the excitation light was found to be low outside the hemispherical portion of the simulated phantom. As a consequence, simulated heterogeneity present in the cylindrical portions did not have an effect on the tomographic reconstruction. Therefore, heterogeneities were simulated only in the hemispherical portion of the breast-shaped geometry. Therefore, in the breast geometry, Ω corresponds to the hemispherical portion and V(Ω) corresponds to its volume (see Eq. 8). In order to position the lumps, a random number between 0 and +5 assigned the z-position of the center of the lump, and two random numbers between -5 and +5 assigned its x and y-position. Only if the resulting center of the lump (x, y, z) resided within the hemisphere, was the heterogeneity counted as one of N lumps, where the value of N was fixed to 100.

Lumps were generated for all endogenous (μ_axi, μ_ami, μ_sx , and μ_sm) and exogenous (μ_axf and μ_amf) optical parameters. Lumps were first independently generated at prescribed percentages for endogenous properties of μ_axi and μ_sx, and the exogenous property of μ_axf. The lump strengths in the optical parameters of μ_ami and μ_amf were taken to be a factor of that in μ_axi and μ_axf respectively. The factor essentially reflects the wavelength-dependent absorption owing to tissue chromophores and can be thought of as the ratio of bulk extinction coefficients of tissue at excitation and emission wavelengths.

Similarly the lump strength in the optical parameter of μ_sm was taken to be a factor of that in μ_sx. The factor effectively represents the ratio of bulk wavelength dependent scattering efficiency at excitation and emission wavelengths. In this study, we simulated indocyanine green (ICG) as the contrast agent. Thus, the multiplication factors were taken considering the absorption and emission spectra of ICG at 780 nm (excitation) and at 830 nm (emission) wavelengths, respectively. The factors for relating lumps in optical properties at the emission wavelength to those independently generated for lumps in optical properties at the excitation wavelength are given in Table 2. For example, in order to generate lumps in the absorption coefficients due to endogenous chromophores, the spatial distribution and strength of lumps in μ_axi were first calculated using Eq. (7), and the strength of these lumps were then multiplied by the corresponding factor from Table 2 to obtain lumps in μ_ami. For endogenous lumpy object models, the forward solution was solved in the presence of lumps in (μ_axi , μ_ami , μ_sx, and μ_sm) as well as in the presence of a prescribed fluorescent target. For exogenous lumpy object models, the forward solution was solved in the presence of lumps in (μ_axi, μ_ami, μ_sx , and μ_sm) and (μ_axf and μ_amf) as well as in the presence of a fluorescent target.

Table 2. Constant factors to obtain lumps at emission wavelength from the generated lumps at excitation wavelength.

View Table | View all tables in this article

The fluorescent target was taken to be a spherical volume of 1 cm³ centered at spatial coordinates (0.5599, -2.4707, 2.4968) with higher contrast than the background.

Generation of the optical property heterogeneity map and solution of the forward solution was conducted on a workstation with Intel Pentium IV 2.80 GHz and 1.0 GB RAM. The calculation of the volume integral in Eq. (8) was a computationally time consuming step and has been performed using a Simpson quadrature method for numerical integration provided in MATLAB^®.

The background tissue optical properties are provided in Table 3. The target was assumed to have the same optical properties as the background except for the values of μ_axf and μ_amf.

The values of μ_axf and μ_amf for the target were 100, 50 and 25 times more than that of the average background values, and are presented as three different case studies in this work. While actual TBR values may be ten-fold or greater, we employ these ratios assuming molecular targeting of agents. Twenty-five point sources and 128 detectors on the hemispherical surface were employed to obtain the synthetic measurements, as shown in Fig. 2.

Table 3. Average background optical properties. Also given are the parameters used in equations (2), (3), and (4). The optical properties used in the simulations are similar to breast tissue optical properties reported in literature (see Table 1).

View Table | View all tables in this article

 

Fig. 2. The positions of laser sources and detectors on the breast phantom. The stars are the positions of the sources and the circles are the positions of the detectors.

Download Full Size | PDF

2.3 Inverse imaging algorithm

The starting point for formulating the reconstruction as an optimization problem begins with defining the error function E(μ_axf),

subject to the constraint {I ≤ μ_axf ≤ u}, where l is the lower and u is the upper bound of μ_axf. For optical tomography problems, a range, i.e., the lower and upper bounds, is specified for the optimization variables based on physical consideration. In above equation cal denotes the value calculated by the forward problem; mea denotes the synthetic measured values and N_B is the number of detectors. The superscript ∗ denotes the complex conjugate of the complex number Z_P. Z_P is comprised of the referenced fluorescent amplitude, I_ACrefp and the referenced phase shift θ_refp measured at boundary point, p, in response to point illumination. Specifically, the referenced measurement at boundary point p can be given by:

The reference scheme was used in the synthetic data to remove the influence of the unknown source strength [28]. In this scheme, the emission fluences, Φ_m, from multiple collection fiber locations were referenced with respect to the excitation fluence data from the same fiber as given in the following equation

Where p is the position of collection on the surface of the phantom. Thus, the relative phase shift and AC ratio were calculated from the difference in phase [i.e., (∆θ_m)_p =(θ_m)_p -(θ_x)_p] and the ratios of amplitude [i.e., (∆I _ACm)_p = (I _ACm)_p /(I_ACx)_p].

The penalty modified barrier function method with simple bounds constrained optimization technique PMBF/CONTN was used for reconstruction and is presented in detail elsewhere [10]. In the PMBF/CONTN method, the modified barrier penalty function, M (termed hereafter as the barrier function), is used to incorporate the constraints directly within the optimization variable:

where η is a penalty/barrier term;f is a logarithmic function; N is the number of nodal points; λ^l and λ^u are vectors of Lagrange multipliers for the bound constraints of the lower and the upper bounds, respectively. From Eq. (12) one can see that the simple bounds are included in the barrier function M.

The PMBF/CONTN method is performed in two stages within an inner and an outer iteration. The outer iteration updates the Lagrange multiplier λ and the barrier parameter η using the formulations presented elsewhere [10]. Using the calculated values of the parameters λ and f, the constrained truncated Newton with trust region method is applied to minimize the penalty barrier function described in Eq. (12). Once satisfactory convergence within the inner iteration is reached or a specified number of inner iterations have been exceeded, the variables describing PMBF convergence and fractional change in the outer iteration are recalculated to check for convergence. If unsatisfactory, the Lagrange multipliers, λ, and the barrier parameter, η, are updated and another constrained optimization is started within the inner iteration. Herein we use the modified method of Breitfeld and Shanno [29] for initializing the Lagrange multipliers without a priori information. This approach effectively removes the need to choose an appropriate initial value of the Lagrange multipliers based upon ground truth, which, in actual patient imaging, may or may not be available. This algorithm is novel in the sense that the constrained have been included within the optimization variable and no a priori information is required for the regularization parameters. The algorithm has been validated from experimental data employing planar and point illumination/collection geometries [10].

The image reconstruction computations were performed on a LINUX workstation with AMD Opteron 250, 2.4 GHz and 4.0 GB RAM. For image reconstruction, synthetic measurements of emission photon density wave amplitude and phase delay were obtained by solving the forward problem in the presence of the fluorescent target and various lumpy optical backgrounds. These synthetic measurements were input into the PMBF/CONTN inverse algorithm to recover the exogenous optical property (μ_axf) distribution for tumor detection task. In this initial study for detecting the target in the reconstructed tissues, we use the optical property map of the reconstructed μ_axf and the a priori knowledge about the location of the actual target. If the reconstructed map shows a single structure whose centroid lies within 1 cm radius of the actual target’s centroid, then we look for any other structures present in the domain. If no other structures are present, then it is deemed that the target was detected successfully. If other structures are present, but they are located outside a 2 cm radius of the centroid of the actual target, then successful detection of the target is accomplished in the presence of artifacts. However, if no structure is present within 1 cm radius of the actual target’s centroid, then the reconstruction algorithm has failed to detect the target. Another case of failure occurs when a structure is reconstructed near the actual target’s location but is less than half of the actual target’s volume. In all of these cases further investigation is recommended to define the criteria for the detection tasks.

2.4 Figures of merit for image analysis

In order to assess the performance of the PMBF/CONTN algorithm for the detection tasks, both qualitative (visual) and quantitative (centroid of the reconstructed target and RMSE) measures were used. The general mesh viewer, GMV, software was used to plot the reconstructed images.

The centroid of the reconstructed target, r _c , was calculated using the following expression:

Where index i represents the nodes present in the reconstructed target. The term N_T denotes total nodes present and r _i denotes the coordinates of the ith node in the reconstructed target. The centroids of the reconstructed targets can be compared with the centroid of the actual target to assess the detection accuracy of the algorithm.

Similarly, the RMSE values were calculated using the expression:

Where N is total number of nodes in the domain. Superscript calc denotes the values obtained using reconstruction algorithms; and actual denotes the actual distribution of μ_axf which was used to generate the synthetic image data sets. RMSE values of the reconstructed target are used to assess the relative accuracy of the algorithm for different lump strengths and TBR.

3. Results and discussion

3.1 Generation of lumpy backgrounds

The generation of lumpy backgrounds required ~40 min of CPU time for each independent optical property parameter (with a predefined lump strength value). Figs. 3(a) through 3(c) illustrate movies of cross sections of an example distribution of μ_axi, μ_sx, and μ_axf for the 25% exogenous and endogenous lumpy background model in the absence of a fluorescent target.

 

Fig. 3. Movies of the lumps in endogenous and exogenous optical properties using Lumpy object model. The lumps in μ_axi (1.08 MB) (a), μ_sx (1.07 MB) (b), and μ_axf (1.10 MB) (c) are shown as the cutplanes to the breast geometry (Fig. 1) parallel to yz-plane. The snapshots shown above are cutplanes passing through x=0.20 cm. The spread of the lumps are 5mm and there are 100 lumps uniformly generated in the hemispherical volume. The lumps have the strength values of 25% of the average background values of μ_axi , μ_sx , and μ_axf given in Table 3.

Download Full Size | PDF

3.2 Tomographic image reconstruction

The “ground truth” target position defined by TBR of 100:1 absorption (μ_axf) in the phantom is shown in Fig. 4.

 

Fig. 4. Actual distribution of exogenous optical property of the target inside the breast phantom. The colorbar shows the values of μ_axf in cm^-1. The target has μ_axf value 100 times more than that of the background. The Fig. shown is a cut plane parallel to yz-plane and passing through point x=0.5. The lumps in the background are not shown here.

Download Full Size | PDF

Fig. 5 illustrates the PMBF/CONTN recovery of absorption coefficient owing to fluorophore from synthetic data with lumpy endogenous background (i.e. lumps in μ_axi , μ_ami , μ_sx , and μ_sm) of 1% lump strength. Reconstructed images for other lump strengths are not shown here for brevity, but the fluorescent target was recovered for lump strength up to 100%. We found that reconstructed images show the insensitivity to the various lump strength values and can be easily reconstructed for 1%-100% lumps in the background. Each image reconstructed with PMBF/CONTN required ~3 hr of CPU time. The centroid and the mean displacement of the reconstructed targets are given in Table 4 and 5, respectively, and show that the reconstructed target locations have little error, indicating good detection capability of the PMBF/CONTN in the presence of endogenous lumps.

The RMSE values as a function of endogenous lump strength for the targets reconstructed are shown in Fig. 6. The figure shows that the changes in RMSE values are small for different lump strengths. However, the small change in RMSE is consistent with the modest change in boundary measurements owing to change of endogenous optical property.^[17,30] It should be noted that the RMSE values for different lump strengths are used to compare the relative error in the image reconstruction task.

 

Fig. 5. (1.06 MB) Movie of the recovered distribution of exogenous optical property in the presence of 1% lumps in endogenous optical properties using the PMBF/CONTN inverse algorithm. The colorbar shows the values of μ_axf in cm^-1. The movie shows cut planes parallel to xz-plane. The black spherical mesh represents the actual target.

Download Full Size | PDF

 

Fig. 6. The root mean square error (RMSE) values of the reconstructed images in the presence of endogenous lumps in the background. The target-to-background ratio (TBR) is 100:1 and the RMSE is calculated with respect to the actual distribution of the fluorophore shown in Fig. 4.

Download Full Size | PDF

Table 4. Centroid of the reconstructed targets with varying lump intensities. The coordinate dimensions are in centimeters.

View Table | View all tables in this article

Table 5. The mean displacement of the reconstructed targets with respect to the actual target’s centroid. The dimension is in centimeters.

View Table | View all tables in this article

In contrast to the case of endogenous optical property heterogeneity, as the lump strength of exogenous contrast becomes comparable to that of the target, the task of discriminating the target from the background becomes more difficult. Except for the case where TBR contrast ratio was 100:1, we could not reconstruct the target in a background containing exogenous and endogenous lump strengths up to 100% as shown for endogenous lumps alone. The target could not be reconstructed for lump strength values exceeding 50% in the case of TBR of 50:1 and 25:1. It should be noted that in case of 50:1 TBR with 100% lumps in the background, we did see a reconstructed volume at the position of actual target, but since its volume was less than that of the actual target’s volume, our criteria did not enable definitive identification.

Figure 7 shows movies of the PMBF/CONTN recovery of absorption coefficient in the presence of 1% endogenous and exogenous lumps in the background for the three cases of TBR of 100:1, 50:1, and 25:1. Figures illustrating recovery at other lump strengths are not shown here for brevity. The centroid and the mean displacement of the reconstructed targets are given in Table 4 and 5, respectively, for all the three cases of TBR values (100:1, 50:1 and 25:1).

The RMSE values were calculated as a function of TBR and are plotted in Fig. 8 for each case of endogenous and exogenous lump strength ranging from 1% to 50%. The general trend is that the RMSE value increases as the contrast decreases for each lump strength. This can be explained by the fact that as the contrast decreases the target has lower fluorophore concentration and thus the background lumps in exogenous optical properties are more pronounced with respect to the target. Thus, the inverse algorithm becomes more prone to fail, which it does at lump strength values exceeding 50% in the case of 50:1 and 25:1 TBR values.

 

Fig. 7. Movies of the recovered distribution of exogenous optical property in the presence of 1% lumps in endogenous as well as exogenous optical properties using the PMBF/CONTN inverse algorithm. Figs. (a) (1.25 MB), (b) (1.26 MB), and (c) (1.27 MB) are respectively for TBR values of 100:1, 50:1, and 25:1. The colorbar shows the values of μ_axf in cm^-1. The movies shown are cut planes parallel to xz-plane. The black spherical mesh represents the actual target.

Download Full Size | PDF

 

Fig. 8. The root mean square error (RMSE) values of the reconstructed images in the presence of endogenous as well as exogenous lumps in the background. The plots for all the three cases of TBR values are shown and the RMSE is calculated with respect to the actual distribution of the fluorophore including the target distribution shown in Fig. 4 and the exogenous lumps in the background. The values for 100% are not shown here since we could not reconstruct the target for the TBR values of 50 and 25.

Download Full Size | PDF

It is important to note that in all of the above reconstruction tasks, the inverse algorithm had no prior knowledge about the lumps. In fact, the inverse algorithm was provided with only the average background values of the optical properties. This is also practically significant because in the real world nothing will be known about the spatial distribution of the optical properties of the tissue to be imaged. Since the approximate average values of the optical properties can be assessed through FDPM spectroscopy measurements, the average background values may be an appropriate input into the inverse imaging algorithm.

Finally, the results contained herein represent single cases of tomographically reconstructed targets in a background of lumps. In order to perform OAIQ for target detection and estimation tasks, random generation of background lumps across hundreds or thousands of cases need to be simulated in order to generate a statistical number of reconstructed images. As a result, we will be able to generate RMSE for lumpy backgrounds with varying lump strengths in order to ascertain error bars in each of these strength value points. This will be helpful in determining the error bars on Fig. 6 in order to better understand its behavior. In this work, we demonstrate the feasibility of generating images from lumpy backgrounds without a priori information. With the ability to reconstruct targets in the presence of lumpy backgrounds, we are now engaged in OAIQ for optical tomography.

4. Conclusion

In fluorescence-enhanced optical tomography, the natural anatomical background is one source of randomness in the task of detecting a molecularly targeting fluorescent agent. Another source will be the lower level expression of disease markers throughout normal tissue which will be responsible for a “background” heterogeneity that also may impact image quality. Computer algorithms for analysis of task performance of fluorescence imaging require a method of simulating natural anatomical background of endogenous optical properties and the heterogeneous background expression of disease markers; as well as a robust image recovery program which does not require the use of a priori information. Herein, we provide a preliminary assessment of the image reconstruction algorithms and lumpy background model as method to establish the feasibility for using OAIQ tools for assessing image performance.