|Home | About | Journals | Submit | Contact Us | Français|
Reconstruction algorithms are presented for two-step solutions of the bioluminescence tomography (BLT) and the fluorescence tomography (FT) problems. In the first step, a continuous wave (cw) diffuse optical tomography (DOT) algorithm is used to reconstruct the tissue optical properties assuming known anatomical information provided by x-ray computed tomography or other methods. Minimization problems are formed based on L1 norm objective functions, where normalized values for the light fluence rates and the corresponding Green’s functions are used. Then an iterative minimization solution shrinks the permissible regions where the sources are allowed by selecting points with higher probability to contribute to the source distribution. Throughout this process the permissible region shrinks from the entire object to just a few points. The optimum reconstructed bioluminescence and fluorescence distributions are chosen to be the results of the iteration corresponding to the permissible region where the objective function has its global minimum This provides efficient BLT and FT reconstruction algorithms without the need for a priori information about the bioluminescence sources or the fluorophore concentration. Multiple small sources and large distributed sources can be reconstructed with good accuracy for the location and the total source power for BLT and the total number of fluorophore molecules for the FT. For non-uniform distributed sources, the size and magnitude become degenerate due to the degrees of freedom available for possible solutions. However, increasing the number of data points by increasing the number of excitation sources can improve the accuracy of reconstruction for non-uniform fluorophore distributions.
Bioluminescence tomography (BLT) and Fluorescence tomography (FT) are imaging methods which provide monitoring of biological processes at a molecular level and can be used for noninvasive visualization of disease progression and response to treatment [1–4]. Optical molecular probes can be used to label cells or proteins and the light they emit can be detected at the surface giving information about the distributions of specific genes or proteins associated with these sources . However, due to multiple light scattering in tissue, the image does not reflect the true distribution of the bioluminescence sources and fluorophore concentration. BLT and FT are inverse problems that attempt to reconstruct the true distribution of the bioluminescence and fluorescence sources from the measured light.
BLT and FT are ill-posed and ill-conditioned problems [6,7]. The limited number of measurements obtained from the boundary compared to the number of unknowns or variables precludes a unique solution to the problem . To reduce the ill-posedness, multi-spectral measurements can increase the measured data for BLT [8–10]. The number of unknowns can be reduced by using a priori information about permissible source regions [11–13] or solving for a rough estimate of an optimal permissible source region . Different researchers have used sparsity regularization [15–17] for reconstructing localized and sparse sources. In our previous work , an iterative minimization solution based on the L1 norm was used to shrink the permissible regions where the sources are allowed by selecting points with higher probability to contribute to the source distribution. This provides an efficient BLT reconstruction algorithm with the ability to determine relative source magnitudes and positions in the presence of noise. However, the reconstruction of the total source power was not accurate due to the use of the regularization penalty term and the non-optimized choice of the permissible region. In addition, it was not possible to use the algorithm for reconstructing large distributed sources as the size of the final permissible region is determined a priori.
In this work, an algorithm to reconstruct distributed as well as sparse bioluminescence sources has been developed. The algorithm, based on diffusion theory [19–26], used the diffuse optical tomography algorithm we have developed previously  to reconstruct the tissue optical properties assuming anatomical information obtained by another imaging modality such as x-ray CT. The reconstruction problems mentioned above have been solved in the current algorithm through the following: First, the objective function of the minimization algorithm avoids the use of a regularization penalty term. Second, the objective function is formed from the summation of the first norm of the difference between the multi-spectral measured light fluence rate and that calculated using the Green’s function at different wavelengths. The light fluence rate and the corresponding Green’s function at each wavelength are normalized. This normalization makes contributions from all wavelengths comparable instead of having the objective function dominated by the longer wavelength terms where the signal is largest. . Third, the permissible region is initially chosen to be the whole object and then iteratively shrunk by scanning all possible combination of points inside the object. The optimized permissible region that gives the best result for the reconstructed source corresponds to the global minimum of the objective function. This provides an objective choice of the best solution and enables the algorithm to reconstruct distributed as well as sparse sources.
Different algorithms have been developed to solve the FT problem. Some groups emphasize the importance of using a priori information [27–29] about the background optical properties of the heterogeneous medium. Other groups used a priori information [30–32] about the fluorophore distribution in different tissues by employing “hard” or “soft” constraints to improve the reconstruction accuracy. Recent simulation studies [33,34] have suggested that both types of a priori information are essential for the accurate recovery of the fluorophore concentration in small inclusion embedded in a heterogeneous medium .
A new algorithm for FT has been developed here to reconstruct distributed and localized fluorophore distributions. The total number of fluorophore molecules is recovered with good accuracy without a priori information about the spatial distribution. First the optical properties of different tissues at the excitation and emission wavelengths are reconstructed using our DOT algorithm . The emission light fluence rate and the corresponding Green’s function for every excitation source are normalized. An algorithm similar to that described above for BLT based on the L1 norm and adaptive shrinking of the permissible region has been used for reconstruction of the fluorophore concentration
The two algorithms for BLT and FT rely on solving the diffusion equation numerically using the method of finite elements (FE). The object used in the simulation represents an idealized cross-section of a mouse abdomen. To satisfy some realistic conditions for in vivo imaging, the illumination of the object for DOT and FT is done using point sources embedded in the animal support to illuminate the posterior surface of the animal (i.e. 90 degrees of the circumference) while transmission and emission measurements are made from the anterior and lateral surfaces (i.e. 270 degrees of the circumference).
where is the light fluence rate at position and wavelength λ. The isotropic source distribution is given by. The spatial distribution of the tissue optical properties at wavelength λ is given by the absorption coefficient and the diffusion coefficient, where the diffusion coefficient is defined by
where is the reduced scattering coefficient , is the scattering coefficient and g is the anisotropy factor. At the boundary, the light fluence rate satisfies the equation
where is a unit vector pointed outwardly normal to the surface, and A is derived from Fresnel’s law as 
where the critical angle and and n is the tissue refractive index. Equation (1) can be solved numerically using the method of finite elements (FE), and the discrete version of Eq. (1), assuming a total number of nodes N in the mesh, is given by 
where is the FE matrix of the diffusion operator on the left hand side, and is the FE assembly matrix multiplied by the source resulting from the FE discretization. The light fluence rate is obtained from Eq. (5) and is given by
where the Green’s function of the system gives the light fluence rate at any point due to a unit excitation source at any other point and is defined as
We assume that the light fluence rate at the boundary can be derived from the emission images captured by the CCD camera  and these data will be used for the source reconstruction. The light fluence rate at the detector positions on the boundary is given by
where is the light fluence rate at the detectors and is a matrix that gives the values of the Green’s function at the detector positions. The detector index is d and if the detector points do not coincide with the edge points of the FE mesh, interpolations are used to get the values of the Green’s function at the detector positions. The estimation of the BL source from Eq. (8) is an ill-posed problem because the number of data points is less than the number of unknown source values. To increase the number of data points and reduce the ill-posedness, bioluminescence images can be acquired at different wavelengths. The number of variables can be reduced by restricting the source to a permissible region. In our previous work , an algorithm for BL reconstruction was developed that used multi spectral measurement and iterative update to the permissible region to reconstruct sparse sources. The algorithm minimized an objective function which was a superposition of the first norm of the difference between the calculated and measured boundary light fluence rate and the zero norm of the source. The number of source points in the final iteration was pre-determined so that only few sparse sources could be reconstructed. A penalty term was essential in the iterative solution to guide the minimizer to choose points deeper inside the object instead of points close to the boundary. The weight of the penalty term will affect the accuracy of the reconstructed source location, strength and total power.
Here we propose an improved algorithm for BL reconstruction. Like our earlier method it relies on DOT to reconstruct the tissue optical properties making use of the anatomical information available from x-ray CT or other techniques . These optical properties are used to calculate the Green’s functions described above. The first improvement is the use of normalization: the light fluence rate and the Green’s function at each wavelength are divided by the maximum value of the light fluence rate at that wavelength. Without normalization, the longer wavelength terms will dominate the objective function due to the smaller values of the absorption coefficient and the minimizer will be guided mainly by these data. Normalization effectively increases the efficiency of the collected data points by making the contributions from all wavelengths equal and so more useful for properly guiding the search direction and permits the second improvement - reconstruction without a penalty term. The third improvement is the selection of the optimum source distribution that minimizes the objective function. As in our previous method, for every iteration, the source points are sorted in order of descending power and a fixed fraction of the points is deleted. Unlike our previous algorithm where this process was continued until the number of source points was equal to the number of data points, the process now continues until there are only a few points. The objective function is calculated at each iteration and the source distribution that minimizes the objective function is taken as the best solution. The minimization problem is given in Eq. (9), while the description of the algorithm is given in Table 1 .
In this paper we extend the ideas described above and in  to the problem of fluorescence tomography. Two coupled diffusion equations are used to describe the propagation of excitation and emission light in tissue:
where and are the excitation and emission light fluence rates, respectively, and are the diffusion and absorption coefficients at the excitation wavelength, and are the diffusion and absorption coefficient at the emission wavelength, is the absorption coefficient of the fluorophore which is proportional to its concentration , , where ε is the extinction coefficient of the fluorophore, C is the concentration, and η is the fluorescence quantum yield. As a specific example, we consider the fluorophore ICG which has excitation and emission maxima at 785 and 830 nm. The extinction coefficients at these wavelengths are 130,000 and 22,000 molar−1mm−1, respectively  and the quantum yield is 0.016 . The boundary condition for the excitation and emission is a Robin-type boundary condition as given in Eq. (3). The starting point for our algorithm is the solution of the DOT problem at both the excitation and emission wavelengths. These estimated optical properties are then used for forward calculations of the coupled diffusion equations.
Equation (10) is solved numerically using the finite element method where its discrete version is given by
where is the FE matrix of the diffusion operator at the excitation wavelength in the left hand side, and is the FE assembly matrix multiplied by the source resulting from the FE discretization. The Green’s function of the system at the emission wavelength for Eq. (11) is given by
The light fluence rate at the emission wavelength is given by
The reconstruction problem is to estimate the concentration of the fluorophore using the measurement of the light fluence rate at the emission wavelength. For each source used for the excitation, the corresponding data for the light fluence rate at both excitation and emission wavelength are measured and then a minimization problem of the difference between the calculated and measured light fluence rate can be formulated as
where in the number of excitation sources, is the light fluence rate at the excitation wavelength due to source i, is the detector index for detectors associated with the source i, and and are the normalized measured light fluence rate and the corresponding normalized Green’s function. The same technique for shrinking the permissible region described above for BLT is also used for FT. In each iteration a fraction of points with the lowest fluorophore concentration is removed. The description of the FT algorithm is given in Table 2 .
The 2D object (Fig. 1 ) used in the simulations is 25 mm in diameter and represents an idealized cross-section of a mouse abdomen where skin (E1), bowel (E2), kidneys (E3, E4), spinal column (E5), and adipose tissue (E6) are identified and assigned realistic wavelength-dependent optical properties. It is a simplified version of Fig. 6(a) in Ref. . It is assumed that different regions are distinguished through CT x-ray imaging or other imaging technique and all elements within each region have the same optical properties. In Fig. 1, the source-detectors setup used in the simulation is shown. Five point sources embedded in the animal support table illuminate the posterior surface of the animal (i.e. 90 degrees of the circumference) while transmission measurements for DOT are made from the anterior and lateral surfaces (i.e. 270 degrees of the circumference) using 121 detectors corresponding to pixels on a CCD camera. Emission from the anterior and lateral surfaces is used as input for, BLT, and FT and the same source locations are used for fluorescence excitation. This geometry allows all imaging data including x-ray CT to be acquired without moving the animal. The total number of detector readings Nd = 5*121 = 605 for DOT and FT at excitation wavelength 785 nm and emission wavelength 830 nm. The total number of detector readings used for the BLT is Nd = 3*121 = 363 where the wavelengths used are 580 nm, 600 nm, and 620 nm.
In the first step, the optical properties of all tissues at the desired wavelengths are reconstructed using the DOT algorithm described in . It is assumed that an accurate segmented image is provided by x-ray CT, MR or by a deformable atlas of mouse anatomy. All elements within each region in the object are assumed to have the same optical properties. The light fluence rate at the boundary is written as a Taylor series expansion around an initial guess corresponding to a homogenous object with the same optical properties in all regions. The fluence rate is approximated by the first three terms including the first and second order derivatives which are calculated by the direct method and give the change of light fluence rate at the boundary due to small changes in the tissue optical properties in each region. The first and second order derivatives are used in an iterative algorithm to reconstruct the tissue optical properties. For this simulation forward finite element calculations of the light fluence rate were performed using the true optical properties and two per cent random Gaussian noise was added to obtain the data for reconstruction. The DOT algorithm recalculates the first order derivatives (Jacobian) around the new initial guess of the scattering and absorption coefficients obtained by the last iteration. The recalculation of the Jacobian using the up-to-date values of the optical properties is more efficient than using only one Jacobian as it modifies the search direction and accounts for nonlinearities of the dependence of the light fluence rate on tissue optical properties. In Table 3 , the true and reconstructed scattering coefficients and, respectively are presented while Table 4 shows the true and reconstructed absorption coefficients and, respectively. The true scattering and absorption coefficients were calculated using the data reported in [42,43]. The maximum relative error obtained for scattering coefficients is 15% at 580 nm while the maximum relative error obtained for the absorption coefficients is 18% at 785 nm. The total number of iterations used is 40.
The true optical properties and the true source distribution were used to calculate the bioluminescence fluence rate and 2% Gaussian noise was added to simulate realistic data. The BLT reconstruction algorithm used the optical properties in Tables 3 and and44 recovered by the DOT algorithm. The matrix elements corresponding to the 121 detectors are used for both the Green’s functions and the simulated light fluence rates.
The objective function is minimized using the fmincon function of Matlab and the source is constrained between 0 and 1.5 nW/ mm2. The results are not sensitive to the value of the upper constraint but selecting a biologically plausible limit reduces the range of possible solutions. The first and second derivatives of the objective function are provided to fmincon to enhance speed. The first derivative used in the model for the objective function is given by, where is a column vector of size equal to the number of detectors and its elements equal 1 for , −1 for , and 0 for where the first norm is not differentiable. The second order derivative used is zero.
The permissible region is initially chosen to be the whole object excluding the points corresponding to one transport length from the surface. At every iteration, the permissible region is slightly shrunk by an arbitrary factor β. For a permissible region initially containing 1389 mesh points that is finally shrunk to 10 points in 60 iterations, β = (1389/10) ^ (1/60) = 1.0857.
Figure 2 shows the value of the objective function obtained in every iteration as a function of the number of source points in the permissible region for the true source distribution illustrated in Fig. 3(a) . It is dimensionless since it is a summation of the first norm of normalized terms. The best solution to the reconstructed BL source is obtained when the objective function records its minimum value 0.3385 when the permissible region contains 95 Points at iteration 33.
Figure 3 shows the results of the BLT reconstruction algorithm for different source distributions. The source distribution in 2D space is given in nW/mm2 which implies that integration over z direction has been performed which is equivalent to multiplying by a unit length of 1 mm. The total source power is then obtained by integration over object area. The actual distribution in Fig. 3(a) is two 3mm diameter sources of equal and uniform magnitude 1 nW/mm2 located in the right and left kidneys. The total power of the left sources is 11.06 nW and the right is 10.59 nW. The total power of the reconstructed left and right sources in Fig. 3(d) are 11.11 and 10.57 nW, both within 5% of the true value, and the magnitudes of the reconstructed sources range between 0.78 and 1.22 nW/mm2. The accuracy of reconstructed source power is important since it reflects the total number of bioluminescent molecules. The total powers of the actual and reconstructed sources are calculated numerically by integrating the source magnitude over the source area which does not appear as a perfect circle in the figure due to discretization of the mesh. Figures 3(b) and 3(e) show the actual and reconstructed source distributions when a 3 mm source of magnitude 0.5 nW/mm2 (total power 5.53 nW) is in the left kidney and a 3 mm diameter source of magnitude 1 nW/mm2 (total power 10.59 nW) is in the right. The total powers of the reconstructed left and right sources are 4.33 nW and 11.84 nW, and their magnitudes range between 0.8 to 1.12 nW/mm2 for the left source and 0.5 to 1.08 nW/mm2 for the right source. Figures 3(c) and 3(f) show the actual and reconstructed source distributions when both sources have magnitude 1 nW/mm2 but the left source is 2 mm in diameter (total power 5.66 nW) and the right is 3 mm (total power 10.59 nW). The total powers of the reconstructed left and right sources are 5.19 nW and 11.46 nW, and their magnitudes range between 0.9 to 1.13 nW/mm2 for the left source and 0.8 to 1.1 nW/mm2 for the right source. The total powers of the actual and reconstructed sources are within 20% in all cases but the magnitudes and sizes of the sources are not reconstructed properly in case Fig. 3(e). The reason is that the magnitude and size of the source are degenerate to the algorithm: as long as the total power is the same, a smaller source with larger magnitude cannot be distinguished from a larger source with smaller magnitude. The shrinking of the permissible region is based on removing points corresponding to small source magnitudes. Therefore, if we have two sources with the same size and different magnitudes, more points will be removed in the vicinity of the source with lower magnitude because it has smaller total power. When the lower magnitude source becomes smaller than the higher magnitude source, the algorithm needs to increase the magnitude of the smaller source above the true value to reconstruct the total power and reduce the objective function.
Figure 4 shows the result of the reconstruction of a large uniform source; the actual source in Fig. 4(a) is of uniform magnitude 1 nW/mm2 and fills the gut region. The total power of the reconstructed sources in Fig. 4(c) is 143 nW which is within 2% of the true value, and the magnitude of the reconstructed sources range between 0.8 and 1.2 nW/mm2. The algorithm can reconstruct distributed sources as well as small sparse sources because the number of points in the permissible region is not specified a priori. The accuracy of the reconstruction is high as long as number of points in the true source is not much greater than the number of measured data. The number of data points for the three wavelengths is 363 and the number of points in the optimized permissible region corresponding to the reconstructed source in Fig. 4(c) is 425 obtained at the iteration 14. However, for a non-uniform source distribution such as that shown in Fig. 4(b) where a 4 mm diameter hot spot with magnitude 1.25 nW/mm2 is located in a uniform background of 0.25 nW/mm2, the reconstruction fails to recover the details of the distribution, although the total power of the source is recovered within 2%. The reason for this failure is, as explained above, the degeneracy of magnitude and size of sources with the same total power.
For these simulations the true optical properties are used to calculate the light fluence rate at the excitation wavelength (785 nm) which is then combined with the true fluorophore concentration, the extinction coefficient, and the quantum yield to calculate the emission light fluence rate at 830 nm using Eq. (14). 2% Gaussian noise is added to the simulated data to obtain the measured light fluence rate used in the minimization problem in Eq. (15). The optical properties recovered by DOT are used to calculate the light fluence rate at the excitation wavelength and the Green’s function at the emission wavelength to yield the emission light fluence rate in the objective function in Eq. (15). The matrix elements corresponding to the 121 detectors are chosen for both the Green’s functions and the simulated light fluence rates for each excitation source, giving 605 readings in total
The minimization problem in Eq. (15) is formulated using the normalized light fluence rate and Green’s functions. The role of this normalization is to adjust the weight of the readings corresponding to each excitation sources so they make comparable contributions to the objective function. The objective function is minimized using the fmincon function of Matlab and the fluorophore concentration is constrained between 0 and 1.5 µmol/ mm2. As for BLT, the choice of this maximum value is not critical as long as it is larger than biologically plausible. The technique described above for BLT for shrinking the permissible region and choosing the best solution corresponding to the optimized permissible region when the objective function has its global minimum is applied.
Figure 5 shows the results of the FT reconstruction algorithm for different fluorophore distributions. The actual distribution in Fig. 5(a) is two 3mm diameter regions of equal and uniform concentration 1 µmol/mm2 located in the right and left kidneys. The total number of molecules in the left region is 11.06 µmol and there are 10.59 µmol in the right. The total number in the reconstructed left and right sources in Fig. 5(d) is 12.28 and 11.12 µmol, within 10% of the true values. The concentration of the reconstructed sources ranges between 1 and 1.16 µmol/mm2 for the left region and 0.8 and 1.2 µmol/mm2 for the right. Figures 5(b) and 5(e) shows the actual and reconstructed fluorophore concentrations when a 3 mm source of uniform concentration 0.5 µmol/mm2 (total number of molecules 5.53 µmol) is in the left kidney and a 3 mm diameter source of uniform concentration 1 µmol/mm2 (total number of molecules 10.59 µmol) is in the right. The total numbers of molecules in the reconstructed left and right regions is 5.74 µmol and 11.94 µmol, and their concentrations ranging between 1.15 and 1.24 µmol/mm2 for the left source and 0.8 and 1.15 µmol/mm2 for the right source. Figures 5(c) and 5(f) shows the actual and reconstructed fluorophore concentrations when both sources have uniform concentration 1 µmol/mm2, but the left source is 2 mm in diameter (total numbers of molecules 5.66 µmol) and the right is 3 mm (total numbers of molecules 10.59 µmol). The total numbers of molecules in the reconstructed left and right sources is 5.73 µmol and 11.72 µmol, and their concentrations range between 1 and 1.15 µmol/mm2 for the left source and 0.8 and 1.06 µmol/mm2 for the right. As in the case of BLT, the total number of fluorophore molecules in the actual and reconstructed sources are within 15% in all cases but the concentrations and sizes of the sources are not reconstructed properly in case Fig. 5(e) when we have multiple regions with different concentrations. The reason is the degeneracy or the degrees of freedom that the current algorithm has for choosing different sizes and concentrations that give the same total number of molecules, as explained above for BLT.
Figure 6 shows the result of the reconstruction of a large uniform fluorophore source; the actual source in Fig. 6(a) is of uniform concentration 1 µmol/mm2 and fills the gut region. The total number of molecules of the reconstructed sources in Fig. 6(c) is 155 µmol which is within 7% of the true value, and the concentration of the reconstructed sources range between 0.8 and 1.2 µmol/mm2. The algorithm can reconstruct distributed sources as well as small sparse sources because the number of points in the permissible region is not specified a priori. The accuracy of the reconstruction is high as long as the number of points in the true source is not much greater than the number of measured data. The number of data points for the 5 excitation sources is 605 and the number of points in the optimized permissible region corresponding to the reconstructed source in Fig. 6(c) is 434. However, for a non-uniform source distribution such as that shown in Fig. 6(b) where a 4 mm diameter hot spot with concentration 1.25 µmol/mm2 is located in a uniform background of 0.25 µmol/mm2, the reconstruction fails to recover the details of the distribution, although the total number of molecules of the fluorophore is recovered within 6%. The reason, as explained above, is the degeneracy of magnitude and size of fluorophore distributions with the same total number of molecules.
Figure 7 show the reconstruction of the non-uniform fluorophore concentration when 16 sources uniformly distributed around the object are used and, for each source, 121 detectors are uniformly distributed in a 270 degree field of view. The total number of data points is 16*121 = 1936. The reconstructed concentration in Fig. 7(b) shows non-uniformity as there is a hot spot with concentration 1 µmol/mm2 which is closer to the actual value of 1.25 µmol/mm2 in Fig. 7(a). Therefore, increasing the number of data points and the field of view can reduce the degeneracy in concentration and size of the fluorophore distribution and improve the reconstruction of non-uniform distributions.
Two numerical algorithms based on finite element methods and the diffusion equation have been developed for BLT and FT. The anatomical structure of the object is assumed to be known and a DOT algorithm is used to reconstruct the tissue optical properties at different wavelengths. The forward model applied with the true tissue optical properties and 2% Gaussian noise has been used to generate realistic data for the reconstruction algorithm for DOT, BLT and FT. The Green’s functions of the object have been calculated using the reconstructed optical properties. The normalized data for the measured light fluence rate and the Green’s function have been used to formulate the objective function for the minimization algorithm. An adaptive shrinking of the permissible region is performed where BL points with small values or FL points with low concentration are removed. The optimum reconstructed BL and FL distributions are chosen to correspond to the iteration with the global minimum of the objective function For BLT the results show good agreement for the total power and location of the actual and reconstructed sources in all cases including a large distributed source. For non-uniform distributions, the location and total power of the sources are reconstructed with good accuracy but the size and source magnitude could not be reconstructed well due to the degeneracy of possible solutions. The FT results show good agreement for the location and the total number of molecules in all cases. For non-uniform concentration, the size of the fluorophore region and concentration could not be reconstructed well due to degeneracy. However, increasing the number of data points can improve the reconstruction of non-uniform distributions.
This work was supported in part by grants from the Ontario Ministry of Research and Innovation postdoctoral fellowship program.