The 2D object (
) 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. [

41]. 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 , 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

*N*_{d} = 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

*N*_{d} = 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 [

18]. 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
, the true and reconstructed scattering coefficients

and

, respectively are presented while
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.

| **Table 3**Scattering coefficients for each region in the heterogeneous object calculated at five different wavelengths. The values of the scattering coefficients are given in (mm^{−1}). |

| **Table 4**Absorption coefficients for each region in the heterogeneous object calculated at five different wavelengths. The values of the absorption coefficients are given in (mm^{−1}). |

4.1 BLT reconstruction

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 and 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/ mm

^{2}. 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.

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
. 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.

shows the results of the BLT reconstruction algorithm for different source distributions. The source distribution in 2D space is given in nW/mm^{2} 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 is two 3mm diameter sources of equal and uniform magnitude 1 nW/mm^{2} 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 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/mm^{2}. 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. and show the actual and reconstructed source distributions when a 3 mm source of magnitude 0.5 nW/mm^{2} (total power 5.53 nW) is in the left kidney and a 3 mm diameter source of magnitude 1 nW/mm^{2} (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/mm^{2} for the left source and 0.5 to 1.08 nW/mm^{2} for the right source. and show the actual and reconstructed source distributions when both sources have magnitude 1 nW/mm^{2} 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/mm^{2} for the left source and 0.8 to 1.1 nW/mm^{2} 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 . 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.

shows the result of the reconstruction of a large uniform source; the actual source in is of uniform magnitude 1 nW/mm^{2} and fills the gut region. The total power of the reconstructed sources in 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/mm^{2}. 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 is 425 obtained at the iteration 14. However, for a non-uniform source distribution such as that shown in where a 4 mm diameter hot spot with magnitude 1.25 nW/mm^{2} is located in a uniform background of 0.25 nW/mm^{2}, 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.

4.2 FT reconstruction

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/ mm

^{2}. 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.

shows the results of the FT reconstruction algorithm for different fluorophore distributions. The actual distribution in is two 3mm diameter regions of equal and uniform concentration 1 µmol/mm^{2} 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 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/mm^{2} for the left region and 0.8 and 1.2 µmol/mm^{2} for the right. and shows the actual and reconstructed fluorophore concentrations when a 3 mm source of uniform concentration 0.5 µmol/mm^{2} (total number of molecules 5.53 µmol) is in the left kidney and a 3 mm diameter source of uniform concentration 1 µmol/mm^{2} (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/mm^{2} for the left source and 0.8 and 1.15 µmol/mm^{2} for the right source. and shows the actual and reconstructed fluorophore concentrations when both sources have uniform concentration 1 µmol/mm^{2}, 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/mm^{2} for the left source and 0.8 and 1.06 µmol/mm^{2} 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 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.

shows the result of the reconstruction of a large uniform fluorophore source; the actual source in is of uniform concentration 1 µmol/mm^{2} and fills the gut region. The total number of molecules of the reconstructed sources in 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/mm^{2}. 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 is 434. However, for a non-uniform source distribution such as that shown in where a 4 mm diameter hot spot with concentration 1.25 µmol/mm^{2} is located in a uniform background of 0.25 µmol/mm^{2}, 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.

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 shows non-uniformity as there is a hot spot with concentration 1 µmol/mm^{2} which is closer to the actual value of 1.25 µmol/mm^{2} in . 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.