## 1.

## Introduction

Diffuse optical tomography in fluorescent mode (fDOT), also called fluorescent molecular tomography when fluorescent probes are used for molecular imaging, is becoming an important pre-clinical non-invasive method for small animal imaging.^{1}2.3.^{–}^{4} fDOT reconstruction algorithms model the propagation of light through the subject under study. The model is usually described by the diffusion approximation, since biological tissues present high scatter for wavelengths in the visible and near infrared spectrum ranges. Under this approximation, the optical properties of tissues are described statistically by the scattering coefficient ${\mu}_{s}$ and the absorption coefficient ${\mu}_{a}$. While analytical solutions exist for media with homogeneous optical properties and simple geometries, more complex numerical solutions are required for media with heterogeneous optical properties or more general spatial domains.^{5}^{,}^{6}

The existence of different scattering and absorption coefficients for different tissue types is a realistic assumption when facing the problem of modeling light transport in living subjects. One solution is to estimate these values using DOT^{7}^{,}^{8} and then perform fDOT.^{6}^{,}^{9}^{,}^{10} However, one cannot simultaneously recover absorption and scattering coefficients in constant-wave mode.^{11} Additionally, prior knowledge of the optical parameters in the frequency domain is necessary for accurate fDOT quantification.^{9} Some authors neglect the different values of the scattering parameter and only estimate the absorption coefficient with DOT. This coefficient is then used to reconstruct the fluorophore distribution.^{12}^{,}^{13} Another approach involves combining previously published values for tissue-specific optical parameters with high-resolution anatomical information provided by CT or MRI; however, there is no consensus in the literature on the actual values of tissue-specific optical parameters, and considerable inter-subject variability is reported.^{14}

An alternative approach^{15} simplifies the problem by assuming homogeneous optical properties and normalizing fluorescence photon density to excitation photon density following the so-called normalized Born approximation. This approach is an adaptation of previous DOT algorithms^{16} and facilitates the use of fDOT for small animal imaging, thanks to its low sensitivity to inhomogeneities and straightforward implementation.^{2}^{,}^{17}

However, only a few studies have assessed the quantification accuracy of the method in the presence of different levels of absorption^{18}^{,}^{19} and, to our knowledge, no studies using either simulations or experiments have assessed the effect of scattering heterogeneities on quantification accuracy. In addition, the effects of tissue-specific optical properties reported in previous DOT-fDOT studies without data normalization cannot be extrapolated to normalized data.^{6}^{,}^{9}^{,}^{10}

The goal of our work was to assess the effect of absorption and scattering coefficients on quantification accuracy. In order to provide a theoretical framework, we obtained analytical expressions to explain the dependence of normalized data on changes in background absorption and scattering for homogeneous media using Green’s function for an infinite space and an infinite slab (Sec. 2.2). In principle, if normalized data are not affected by inaccuracies in optical property modeling, they should be invariant to changes in the scattering or absorption of the medium. We then studied the quantification error using four computer-simulated phantoms based on finite element methods. The first phantom was a homogeneous slab with different background absorption and scattering coefficients (Sec. 2.6.1) that was used to validate the previously obtained theoretical result. The second was a heterogeneous slab with a fluorescent cylindrical region in which absorption and scattering coefficients were defined as two, four and six times the background value (Sec. 2.6.2) to cover the wide range of absorption and scattering coefficients reported in the literature.^{14}^{,}^{20} The third was a computer-simulated mouse model that included different scattering coefficients for lung and liver and was used to assess the relationship between quantification accuracy and the tissue-specific values of the scattering coefficient (Sec. 2.6.3). Image reconstructions were performed with the algebraic reconstruction technique based on two different models: 1) assuming homogeneous media and 2) realistically modeling the heterogeneous optical properties using a Jacobian matrix for which heterogeneous optical properties can be explicitly modeled by the finite element method. The resulting reconstructed images were compared quantitatively.

We hypothesized that the assumption of homogeneous optical properties for modeling diffuse media with a highly heterogeneous scattering coefficient could lead to noticeable quantification errors despite data normalization. Accordingly, we also assessed whether a reconstruction method controlling for heterogeneous scattering coefficients could lead to an improvement in quantification accuracy.

## 2.

## Methods

## 2.1.

### Forward Problem

## 2.1.1.

#### Diffusion approximation

Let ${\mu}_{a}(\mathbf{r},\lambda )$ be the absorption, ${\mu}_{s}^{\prime}(\mathbf{r},\lambda )$ the reduced scattering, and $\kappa (\mathbf{r},\lambda )={\{3[{\mu}_{a}(\mathbf{r},\lambda )+{\mu}_{s}^{\prime}(\mathbf{r},\lambda )]\}}^{-1}$ the diffusion coefficient for a wavelength $\lambda $ at position $\mathbf{r}$ in a domain $\mathrm{\Omega}$. In a diffuse range, which is valid for highly scattering (${\mu}_{s}^{\prime}\gg {\mu}_{a}$) and weakly anisotropic media, photon density is approximated by the diffusion equation^{5}. Hence, in fDOT, the excitation photon density ${\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r},{\lambda}_{\mathrm{ex}})$ at the excitation wavelength ${\lambda}_{\mathrm{ex}}$ and emission photon density ${\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r},{\lambda}_{\mathrm{em}})$ at the emission wavelength ${\lambda}_{\mathrm{em}}$ are given by the solution of a pair of coupled diffusion equations.^{6}^{,}^{21} The excitation photon density is emitted by an external source ${q}_{0}({\mathbf{r}}_{s})$ at a location ${\mathbf{r}}_{s}\in \mathrm{\Omega}$, and the emission photon intensity comes from a fluorescent region characterized by a fluorescent yield $F({\mathbf{r}}_{\mathrm{fl}})$, which accounts for its quantum efficiency, absorption parameter and the concentration of fluorophore. Assuming a time-stationary problem (frequency $w=0$) and that the absorption parameter is not affected by the presence of the fluorophore, excitation and emission photon densities are given by

## (1)

$$-\nabla \xb7\kappa (\mathbf{r},{\lambda}_{\mathrm{ex}})\nabla {\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r},{\lambda}_{\mathrm{ex}})+{\mu}_{a}(\mathbf{r},{\lambda}_{\mathrm{ex}}){\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r},{\lambda}_{\mathrm{ex}})={q}_{0}({\mathbf{r}}_{s})$$## (2)

$$-\nabla \xb7\kappa (\mathbf{r},{\lambda}_{\mathrm{em}})\nabla {\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r},{\lambda}_{\mathrm{em}})+{\mu}_{a}(\mathbf{r},{\lambda}_{\mathrm{em}}){\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r},{\lambda}_{\mathrm{em}})=F({\mathbf{r}}_{\mathrm{fl}}){\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r},{\lambda}_{\mathrm{ex}})\mathrm{.}$$In Eqs. (1) and (2), optical coefficients are different at excitation and emission wavelengths, and should be modeled differently in practice. However, since our main goal was to assess the impact of scattering, we believe that including additional problems such as optical properties mismatch at emission and excitation wavelengths would introduce additional variations. Hence, we opted for this simplified expression

## (3)

$$-\nabla \xb7\kappa (\mathbf{r})\nabla {\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r})+{\mu}_{a}(\mathbf{r}){\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r})={q}_{0}({\mathbf{r}}_{s})$$## (4)

$$-\nabla \xb7\kappa (\mathbf{r})\nabla {\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r})+{\mu}_{a}(\mathbf{r}){\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r})=F({\mathbf{r}}_{\mathrm{fl}}){\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r}).$$^{22}and by applying numerical methods, such as the finite element method, for heterogeneous media and/or general geometries.

^{5}

^{,}

^{23}

From a general perspective, let us define a Green function $\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})$ that solves the heterogeneous problem

## (5)

$$[-\nabla \xb7\kappa (\mathbf{r})\nabla +{\mu}_{a}(\mathbf{r})]\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime}).$$## (6)

$${\mathrm{\Phi}}_{\mathrm{ex}}(\mathbf{r})=\int {\mathrm{d}\mathbf{r}}^{\prime}\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime}){q}_{0}({\mathbf{r}}^{\prime})$$## (7)

$${\mathrm{\Phi}}_{\mathrm{em}}(\mathbf{r})=\int {\mathrm{d}\mathbf{r}}^{\prime}\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})F({\mathbf{r}}^{\prime}){\mathrm{\Phi}}_{\mathrm{ex}}({\mathbf{r}}^{\prime}).$$^{24}

^{,}

^{25}In the finite element implementation, photon density and concentration of fluorophore were approximated by piecewise-linear nodal basis functions. Sources [RHS on Eq. (3)] were modeled as isotropic point sources (located at a depth $1/{\mu}_{s}^{\prime}$ below the surface), which resemble a collimated laser, as described in Ref. 25. Measurements were modeled by a Gaussian kernel centered at the detector location and computed as a linear operator $\mathcal{M}$ acting on the photon density at the boundary of the domain. Thus, measured excitation and emission photon densities at the detector position ${\mathbf{r}}_{d}$ are given by

## 2.1.2.

#### Normalized data: general case

Normalized data for fDOT, that is, the incorrectly termed Born ratio, are given^{15}^{,}^{16}^{,}^{18} by

## (9)

$${U}_{b}({\mathbf{r}}_{d})=\frac{{\mathrm{\Phi}}_{\mathrm{em}}^{\text{meas}}({\mathbf{r}}_{d})}{{\mathrm{\Phi}}_{\mathrm{ex}}^{\text{meas}}({\mathbf{r}}_{d})}=\frac{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})F({\mathbf{r}}^{\prime}){\mathrm{\Phi}}_{\mathrm{ex}}({\mathbf{r}}^{\prime})}{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})q({\mathbf{r}}^{\prime})}\phantom{\rule{0ex}{0ex}}=\frac{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})F({\mathbf{r}}^{\prime})\int \mathrm{d}{\mathbf{r}}^{\prime \prime}\mathbf{G}({\mathbf{r}}^{\prime};{\mathbf{r}}^{\prime \prime})q({\mathbf{r}}^{\prime \prime})}{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})q({\mathbf{r}}^{\prime})},$$## 2.2.

### Normalized Data: Homogeneous Media

In this section, we used Green’s function for homogeneous media to derive analytical expressions that describe the dependence on normalized data with background absorption and scattering coefficients. We theoretically solved two cases of interest: infinite space and infinite slab.

Assuming a point source $q(\mathbf{r})=\delta (\mathbf{r}-{\mathbf{r}}_{s})$ and a fluorescent point $F(\mathbf{r})=\delta (\mathbf{r}-{\mathbf{r}}_{f})$, after substitution in Eq. (9), leads to a simplified form, as follows:

## (10)

$$u({\mathbf{r}}_{d};{\mathbf{r}}_{s},{\mathbf{r}}_{f})\phantom{\rule{0ex}{0ex}}=\frac{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})\delta ({\mathbf{r}}^{\prime}-{\mathbf{r}}_{f})\int \mathrm{d}{\mathbf{r}}^{\prime \prime}\mathbf{G}({\mathbf{r}}^{\prime};{\mathbf{r}}^{\prime \prime})\delta ({\mathbf{r}}^{\prime \prime}-{\mathbf{r}}_{s})}{\int \mathrm{d}{\mathbf{r}}^{\prime}\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}^{\prime})\delta ({\mathbf{r}}^{\prime}-{\mathbf{r}}_{s})}\phantom{\rule{0ex}{0ex}}=\frac{\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}_{f})\mathbf{G}({\mathbf{r}}_{f};{\mathbf{r}}_{s})}{\mathbf{G}({\mathbf{r}}_{d};{\mathbf{r}}_{s})}.$$## 2.2.1.

#### Infinite space

The Green function $\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})$ for an infinite space^{22} is

## (11)

$$\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}\text{\hspace{0.17em}}\mathrm{exp}\text{\hspace{0.17em}}(-\sqrt{3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})}|\mathbf{r}-{\mathbf{r}}^{\prime}|)\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime}|},$$## (12)

$$u({\mathbf{r}}_{d};{\mathbf{r}}_{s},{\mathbf{r}}_{f})=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}\text{\hspace{0.17em}}\mathrm{exp}[-\sqrt{3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})}\phantom{\rule{0ex}{0ex}}(|{\mathbf{r}}_{f}-{\mathbf{r}}_{s}|+|{\mathbf{r}}_{d}-{\mathbf{r}}_{f}|-|{\mathbf{r}}_{d}-{\mathbf{r}}_{s}|)]\phantom{\rule{0ex}{0ex}}\frac{|{\mathbf{r}}_{d}-{\mathbf{r}}_{s}|}{|{\mathbf{r}}_{s}-{\mathbf{r}}_{f}||{\mathbf{r}}_{d}-{\mathbf{r}}_{f}|},$$## (13)

$$u({\mathbf{r}}_{d};{\mathbf{r}}_{s},{\mathbf{r}}_{f})=\frac{3{\mu}_{s}^{\prime}}{2{(2\pi )}^{\frac{3}{2}}}\text{\hspace{0.17em}}\mathrm{exp}[-\sqrt{3{\mu}_{a}{\mu}_{s}^{\prime}}(|{\mathbf{r}}_{f}-{\mathbf{r}}_{s}|+|{\mathbf{r}}_{d}-{\mathbf{r}}_{f}|\phantom{\rule{0ex}{0ex}}-|{\mathbf{r}}_{d}-{\mathbf{r}}_{s}|)]\frac{|{\mathbf{r}}_{d}-{\mathbf{r}}_{s}|}{|{\mathbf{r}}_{s}-{\mathbf{r}}_{f}||{\mathbf{r}}_{d}-{\mathbf{r}}_{f}|}.$$Thus, normalized data are invariant with absorption but almost linearly related to scattering, as the exponential term in Eq. (12) is less dominant (see results, Sec. 3.1).

## 2.2.2.

#### Infinite slab

The Green function $\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime})$ for an infinite slab^{22} is

## (14)

$${{\mathbf{G}}_{\text{slab}}}_{\text{}}(\mathbf{r},{\mathbf{r}}^{\prime})=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}\sum _{n=-\infty}^{\infty}[\frac{\mathrm{exp}(-\sqrt{3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})}{\rho}_{12+{n}^{\prime}})}{{\rho}_{12+{n}^{\prime}}}-\frac{\mathrm{exp}(-\sqrt{3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})}{\rho}_{12-{n}^{\prime}})}{{\rho}_{12-{n}^{\prime}}}]$$## (15)

$$=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}\sum _{n=-\infty}^{\infty}\stackrel{~}{h}({\mu}_{a},{\mu}_{s}^{\prime},\mathbf{r},{\mathbf{r}}^{\prime})\phantom{\rule{0ex}{0ex}}=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}h({\mu}_{a},{\mu}_{s}^{\prime},\mathbf{r},{\mathbf{r}}^{\prime}),$$## (16)

$${\rho}_{12+{n}^{\prime}}^{2}={({x}_{2}-{x}_{1})}^{2}+{({y}_{2}-{y}_{1})}^{2}+{(2\mathrm{nd}+{z}_{2}+{z}_{1})}^{2}$$## (17)

$${\rho}_{12-{n}^{\prime}}^{2}={({x}_{2}-{x}_{1})}^{2}+{({y}_{2}-{y}_{1})}^{2}+{(2\mathrm{nd}+{z}_{2}-{z}_{1})}^{2}$$Thus, the normalized data are given by

## (18)

$$u({\mathbf{r}}_{d};{\mathbf{r}}_{s},{\mathbf{r}}_{f})=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}\frac{h({\mu}_{a},{\mu}_{s}^{\prime},{\mathbf{r}}_{s},{\mathbf{r}}_{f})h({\mu}_{a},{\mu}_{s}^{\prime},{\mathbf{r}}_{f},{\mathbf{r}}_{d})}{h({\mu}_{a},{\mu}_{s}^{\prime},{\mathbf{r}}_{s},{\mathbf{r}}_{d})}\phantom{\rule{0ex}{0ex}}=\frac{1}{{(2\pi )}^{\frac{3}{2}}2\kappa}f({\mu}_{a},{\mu}_{s}^{\prime},{\mathbf{r}}_{s},{\mathbf{r}}_{f},{\mathbf{r}}_{d}),$$## (19)

$$u({\mathbf{r}}_{d};{\mathbf{r}}_{s},{\mathbf{r}}_{f})=\frac{3{\mu}_{s}^{\prime}}{{2(2\pi )}^{\frac{3}{2}}}f({\mu}_{a},{\mu}_{s}^{\prime},{\mathbf{r}}_{s},{\mathbf{r}}_{f},{\mathbf{r}}_{d}).$$Hence, we obtain a relationship similar to that described in the case of infinite space.

## 2.3.

### Jacobian Matrix

The emission photon density ${\mathrm{\Phi}}_{\mathrm{em}}$ is linearly dependent on the fluorophore concentration $F$ under certain assumptions. This can be seen by taking its Fréchet derivative with respect to the concentration of fluorophore and by assuming that a variation of fluorophore $\delta F$ induces a variation of emission photon density $\delta {\mathrm{\Phi}}_{\mathrm{em}}$ while the excitation photon density ${\mathrm{\Phi}}_{\mathrm{ex}}$ (and implicitly the optical parameters) remains invariant and there is no re-emission by the fluorophore,

## (20)

$$\frac{\partial}{\partial \tau}\left[\frac{{\mathrm{\Phi}}_{\mathrm{em}}(F+\tau \delta F)}{{\mathrm{\Phi}}_{\mathrm{ex}}}\right]{|}_{\tau =0}=\nabla {U}_{b}\delta F.$$## (21)

$$\delta {\mathrm{\Phi}}_{\mathrm{em}}({\mathbf{r}}_{d})=\sum _{j}{\int}_{{\mathrm{\Omega}}_{j}}{\mathrm{d}\mathbf{r}}_{j}\mathbf{G}({\mathbf{r}}_{d},{\mathbf{r}}_{j}){\mathrm{\Phi}}_{\mathrm{ex}}({\mathbf{r}}_{j},{\mathbf{r}}_{s})\delta {F}_{j}\phantom{\rule{0ex}{0ex}}=\sum _{j}{\int}_{{\mathrm{\Omega}}_{j}}{\mathrm{d}\mathbf{r}}_{j}\stackrel{~}{\mathrm{\Phi}}({\mathbf{r}}_{j},{\mathbf{r}}_{d}){\mathrm{\Phi}}_{\mathrm{ex}}({\mathbf{r}}_{j},{\mathbf{r}}_{s})\delta {F}_{j},$$^{24}and $\stackrel{~}{\mathrm{\Phi}}({\mathbf{r}}_{j},{\mathbf{r}}_{d})$ is the adjoint field at ${\mathbf{r}}_{j}$ due to a source ${\stackrel{~}{q}}_{0}$ located at the detector position ${\mathbf{r}}_{d}$, which solves an adjoint equation

## (22)

$$[-\nabla \xb7\kappa (\mathbf{r})\nabla +{\mu}_{a}(\mathbf{r})]\stackrel{~}{\mathrm{\Phi}}(\mathbf{r})={\stackrel{~}{q}}_{0}({\mathbf{r}}_{d}).$$## (23)

$${J}_{ij}=\frac{\partial {({U}_{b})}_{i}}{\partial {F}_{j}}=\frac{1}{{\mathrm{\Phi}}_{\mathrm{ex}}^{\text{meas}}({\mathbf{r}}_{d},{\mathbf{r}}_{s})}{\int}_{{\mathrm{\Omega}}_{j}}{\mathrm{d}\mathbf{r}}_{j}\stackrel{~}{\mathrm{\Phi}}({\mathbf{r}}_{j},{\mathbf{r}}_{d}){\mathrm{\Phi}}_{\mathrm{ex}}({\mathbf{r}}_{j},{\mathbf{r}}_{s}).$$^{6}

^{,}

^{21}

^{,}

^{26}

The Jacobian matrix in Eq. (23) was computed using the finite element method, which enables modeling of heterogeneous optical properties.

## 2.4.

### Inverse Problem and Image Analysis

Reconstruction was performed using a randomized algebraic reconstruction technique which has been previously applied in fDOT.^{7}^{,}^{18}^{,}^{27}^{,}^{28} A relaxation parameter of 0.1 and 20 iterations was applied to all the reconstructions.

Image quantification was performed by defining a Region of Interest (RoI) that corresponded to the known target distribution and by taking the average of the voxels within this RoI. The final concentration of fluorophore is expressed as a percentage of the corresponding value in the homogeneous phantom, which was used as a reference. Thus, concentrations become expressed in a relative scale as

where $F$ and ${F}_{\text{homo}}$ are the recovered concentrations of fluorophore corresponding to the heterogeneous and homogeneous model, respectively. The measurement error is also expressed as a percentage, given by $100(F-{F}_{\text{homo}})/{F}_{\text{homo}}$.In addition, profiles were drawn along the reconstructed images as a test of image quality and positional accuracy.

Five percent additive white noise was incorporated into normalized data before reconstruction.

## 2.5.

### Phantoms and Acquisition Protocol

## 2.5.1.

#### Slab phantom

We defined a heterogeneous slab-geometry phantom with two regions, a cylindrical fluorescent region in which optical parameters could be varied and a background region. The slab was 10 mm thick and the cylinder diameter was 5 mm (Fig. 1). A homogeneous phantom was used as a reference.

The acquisition protocol used a configuration of $9\times 9$sources located on the lower plane covering a surface of 12 by $12\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{2}$ (Fig. 1), together with an equivalent configuration of detectors located at the corresponding positions on the upper plane

Computer-simulated data for each experimental phantom were generated using a fine tetrahedral finite element mesh (115 350 nodes and 540 000 tetrahedra). A coarser mesh (55 061 nodes and 253 125 tetrahedra) was used for calculation of the Jacobian matrix, which was resampled into a uniform mesh of $20\times 20\times 20$ for the reconstruction.

## 2.5.2.

#### Mouse model phantom

A heterogeneous model of a mouse that simulates differences in lung and liver scattering was designed based on the Digimouse atlas.^{29} Resembling the experimental set-up proposed in Ref. 1, we simulated a parallel-plate imaging chamber, where the mouse is gently compressed by moving plates filled with intralipid. The model comprised four tissue types: lung, liver, surrounding mouse tissue, and intralipid. It was created by selecting lung and liver regions on the Digimouse model and mapping them into a slab-geometry finite element mesh measuring 12.5 mm in thickness (Fig. 2). The heterogeneous model was created by assigning different ${\mu}_{s}^{\prime}$ for these four regions: $\text{2.12\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for the lung and $0.65\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for the liver as reported in Ref. 14 (for a wavelength of 700 nm); $1.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for the surrounding tissue, which was also used in Ref. 19; and $1.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for the intralipid region selected as a middle value. The absorption parameter was set to ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and kept constant for all models.

The simulated acquisition configuration consisted of $10\times 10$ sources located at the inferior plane covering a surface of $15\times 15\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{2}$, and the same configuration of detectors was used at the corresponding positions on the superior plane.

Finally, two spherical fluorophore spots measuring 2 mm in diameter were simulated and placed inside the right lung and the liver at the same $z$-level ($z=6.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$) (Fig. 2).

## 2.6.

### Quantification Analysis

## 2.6.1.

#### Homogeneous slab phantom

For the homogeneous slab phantom described in Sec. 2.5.1, adopting the parameters ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ as a reference, quantification errors were studied as a function of an $n$-fold change in the background absorption coefficient and an $n$-fold change in the scattering coefficient ($n=0.5$, 1.6, 2.7, 3.8, 4.9, 6). All data sets corresponded to identical spherical fluorophore distributions measuring 2.5 mm in diameter and were reconstructed assuming the reference coefficients.

## 2.6.2.

#### Heterogeneous slab phantom

Heterogeneous slab-geometry phantoms (Fig. 1) were simulated by increasing the scattering coefficient ($2{\mu}_{s}^{\prime}$, $4{\mu}_{s}^{\prime}$ and $6{\mu}_{s}^{\prime}$) or the absorption coefficient ($2{\mu}_{a}$, $4{\mu}_{a}$ and $6{\mu}_{a}$) in the cylindrical region with respect to the background optical parameters (${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$).

All data sets were reconstructed 1) by assuming a homogeneous medium and 2) by modeling the differences in scattering and absorption coefficients. The reconstructed images were quantitatively compared with the reconstruction corresponding to the homogeneous phantom as detailed above.

## 2.6.3.

#### Mouse model phantom

Four reconstructions of this phantom were generated. The first one (reference phantom) made use of the homogeneous model with ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, for both the simulation and the reconstruction. The remaining three were as follows: a) a heterogeneous model reconstructed using the heterogeneous model, b) a homogeneous model with optical parameters chosen as in Ref. 19 (${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$), and c) a homogeneous model with scattering computed as the volume-weighted average of the heterogeneous tissue values (${\mu}_{s}^{\prime}=1.27\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$). The reconstructed images were quantitatively compared with the reconstruction corresponding to the homogeneous phantom, as detailed above.

## 3.

## Results

## 3.1.

### Normalized Data: Homogeneous Media

The dependence of the normalized data on background absorption and scattering was studied for the homogeneous infinite space solution [Eq. (12)]. Normalized data were generated for a point source located at (0,0,0), a point fluorophore located at (0,0,5) and a moving position of the detector along the $x$-axis (located at $(x,0,10)$, $x$ varying from 0 to 6 mm). We found that the exponential term in Eq. (12) is equal to or less than one and decreases slightly when $n$ increases, being equal to one for $x=0$. Hence, the change in normalized data is almost linear with a change in scattering and is almost unaffected by a change in absorption (Fig. 3).

## 3.2.

### Quantification Analysis

## 3.2.1.

#### Homogeneous slab phantom

Quantification errors increased almost linearly with a change in background scattering and increased slightly with a change in background absorption (Fig. 4). For the homogeneous case, the ratio of these two should be linear, as it can be extrapolated from Eq. (18), where the exponential term and the spatial dependence vanished.

## 3.2.2.

#### Heterogeneous slab phantom

Figure 5 shows the relative fluorophore concentration values [Eq. (25)] for changes in absorption and in scattering in the cylindrical heterogeneous region. Reconstructed images are shown in Fig. 6. Quantification errors induced by assuming homogeneous media ranged from $+41$ to $+94\%$ and from $+0.1$ to $-7\%$ for the different values used for scattering coefficients ($2{\mu}_{s}^{\prime}$, $4{\mu}_{s}^{\prime}$, $6{\mu}_{s}^{\prime}$) and absorption coefficients ($2{\mu}_{a}$, $4{\mu}_{a}$, $6{\mu}_{a}$), respectively. Using the reconstruction model that accounted for the previous optical differences led to lower quantification errors in the range $-4$ to $+4\%$, for changes in both scattering and absorption.

The profiles drawn along the $x$-axis on the images in Fig. 6 are also shown (Fig. 7). It seems that positional accuracy is not affected by cylindrical heterogeneity. In contrast, quantification is highly affected by scattering heterogeneity.

## 3.2.3.

#### Mouse model phantom

Quantification results of the relative concentration of fluorophore in the right lung and liver and its ratio are presented in Table 1. In this case, the reference concentration in Eq. (25) is the average between the recovered concentrations in the lung and liver for the homogeneous model. Quantification errors for the lung and liver analyzed separately were in the range $-39$ to $+44\%$ when assuming a homogeneous model. Accounting for the scattering differences improved quantification, with an error of $-7$ to $+7\%$. The quantification error for the ratio between concentrations of fluorophore in the liver and in the lung was $-46\%$ for the homogeneous model, $-13\%$ for the heterogeneous model, and $+6\%$ for the reference model. In the reference approach, a homogeneous model was adopted for both data simulation and reconstruction.

## Table 1

Image quantification of the relative fluorophore concentration [Eq. (25)] in a computer-simulated mouse model with two fluorophore spots, one inside the lung (LU) and one inside the liver (LI). The first row corresponds to simulation and reconstruction using the homogeneous model (homogeneous, μs′=1 mm−1), where the average between the recovered concentrations in LU and LI was used as a reference for the quantification. The second to fourth rows correspond to heterogeneous model simulations reconstructed with the heterogeneous model (second row), the reference homogeneous model (third row), and a homogeneous model with scattering value computed as a weighted average of the scattering values (homogeneous μs′=1.27 mm−1) (fourth row).

Data model | Rec. model | LI | LU | LI/LU |
---|---|---|---|---|

homogeneous | homogeneous | 1,03 | 0,97 | 1,06 |

heterogeneous | heterogeneous | 0,93 | 1,07 | 0,87 |

heterogeneous | homogeneous | 0,78 | 1,44 | 0,54 |

heterogeneous | homogeneous* | 0,61 | 1,14 | 0,54 |

Assuming a homogeneous model led to a decrease in image quality with a distortion from the theoretical spherical shape, as compared with the heterogeneous model (Fig. 8).

Profiles for the reconstructed images in Fig. 8 were drawn along the $x$-axis across the two spherical targets in the lung and liver. As the liver has a smaller scattering coefficient (${\mu}_{s}^{\prime}=0.65\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$) than the background and the lung has a larger one (${\mu}_{s}^{\prime}=2.12\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$), the fluorophore concentration can be both underestimated and overestimated (Fig. 9).

## 4.

## Discussion

Reconstruction based on normalized data has previously been found to be robust in the presence of absorption differences.^{18} To our knowledge, this is the first report presenting evidence of the effect of scattering heterogeneities and quantification accuracy using normalized data. Our results suggest that scattering heterogeneity has a noticeable negative effect on quantification accuracy.

For a homogeneous infinite space, normalized data present a dominant linear relationship with the scattering coefficient while remaining invariant with the absorption coefficient. Accordingly, for a homogeneous slab phantom, a reconstruction assuming incorrect optical parameters leads to quantification errors that increase almost linearly with the scattering coefficient while increasing slightly with the absorption coefficient.

For heterogeneous slab phantoms with different scattering coefficient in a cylindrical fluorescent region, we found that quantification errors induced by assuming a homogeneous medium ranged from $+41$ to $+94\%$, taking quantification from a homogeneous phantom as a reference. Additionally, the presence of regions with absorption heterogeneities led to much smaller quantification errors, ranging from $+0.1$ to $-7\%$, which agrees with the findings of previous reports.^{18} Quantification errors decreased significantly (range, $-4$ to $+4\%$) when using a reconstruction model that accounted for the differences in scattering and absorption coefficients.

For the reconstruction of a computer-simulated mouse model that included different liver and lung scattering coefficient, we obtained quantification errors ranging from $-39$ to $+44\%$ when assuming a homogeneous medium and $-7$ to $+7\%$ when using a reconstruction model that accounted for the differences in scattering. Thus, neglecting the effect of scattering heterogeneities led to a systematic overestimation of the calculated concentration of fluorophore; however, the effect of absorption was remarkably smaller. In more realistic situations, the actual effect may be more complicated.

These results are consistent with those of a previous study that analyzed the effect of modeling errors simultaneously in both scattering and absorption parameters.^{19} The model used in this study was a computer-simulated reconstruction of a mouse torso comprising five regions (liver, heart, lung, bone and background tissues), and relative 2-norm errors only increased 0.74% when a heterogeneous reconstruction model was used. The averaged optical property values in this model were taken from the literature. However using a homogeneous model, errors increased 4.8%. These errors were defined as 2-norm errors that accounted for the entire image and were relative to the error associated with the exact model; consequently, they do not correspond exactly to our quantification errors. Our results correspond to a model that accounted for two organs with distinct scattering (reduced scattering in the range 0.65 to $2.34\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ at 700 nm). Larger errors should be expected when the biological tissue diversity in real animal experiments is taken into account.

From the profiles drawn on the reconstructed images for the heterogeneous slab phantom and for the mouse model phantom, it seems that positional accuracy is not much affected by optical heterogeneity. In contrast, quantification is highly affected by scattering heterogeneity.

The heterogeneous phantoms used in the present study emulate quantification in two organs with different scattering coefficients. Hence, these results cannot be generalized to other arrangements, such as quantification in the presence of organs with a distinct scattering coefficient. Further studies, especially *in vivo* validations, would help to establish the limits of use for normalized data in preclinical studies.

In conclusion, the present study makes use of theoretical results for homogeneous media and computer simulations for heterogeneous models to demonstrate a) that inaccuracies in modeling the scattering parameter may lead to noticeable quantification errors when using a normalized data reconstruction approach and b) that more sophisticated models of the forward problem yield much better results.

## Acknowledgments

This study was supported by Ministerio de Ciencia e Innovación (FPI program, TEC 2008-06715, and CENIT AMIT CEN-20101014), Comunidad de Madrid and European Regional Development Fund ARTEMIS S2009/DPI-1802, and EU-FP7 project FMTXCT-201792.

## References

E. E. Graveset al., “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. 30(5), 901–911 (2003).MPHYA60094-2405http://dx.doi.org/10.1118/1.1568977Google Scholar

V. Ntziachristoset al., “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005).NABIF91087-0156http://dx.doi.org/10.1038/nbt1074Google Scholar

S. Patwardhanet al., “Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice,” Opt. Express 13(7), 2564–2577 (2005).OPEXFF1094-4087http://dx.doi.org/10.1364/OPEX.13.002564Google Scholar

T. LasserV. Ntziachristos, “Optimization of 360 degrees projection fluorescence molecular tomography,” Med. Image Anal. 11(4), 389–399 (2007).1361-8415http://dx.doi.org/10.1016/j.media.2007.04.003Google Scholar

S. R. ArridgeM. HiraokaD. T. Delpy, “A finite element approach for modelling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993).MPHYA60094-2405http://dx.doi.org/10.1118/1.597069Google Scholar

A. Corluet al., “Three-dimensional *in vivo* fluorescence diffusive optical tomography of breast cancer in humans,” Opt. Express 15(11), 6696–6716 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.006696Google Scholar

S. R. ArridgeJ. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Prob. 25(12), 123010 (59pp) (2009).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/25/12/123010Google Scholar

T. Tarvainenet al., “Corrections to linear methods for diffuse optical tomography using approximation error modelling,” Biomed. Opt. Express 1(1), 209–222 (2010).2156-7085http://dx.doi.org/10.1364/BOE.1.000209Google Scholar

Y. Linet al., “Fluorescence diffuse optical tomography with functional and anatomical a priori information: feasibility study,” Phys. Med. Biol. 52(18), 5569–5585 (2007).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/52/18/007Google Scholar

Y. TanH. Jiang, “Diffuse optical tomography guided quantitative fluorescence molecular tomography,” Appl. Opt. 47(12), 2011–2016 (2008).APOPAI0003-6935http://dx.doi.org/10.1364/AO.47.002011Google Scholar

S. R. ArridgeW. R. B. Lionheart, “Non-uniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.23.000882Google Scholar

A. Garofalakiset al., “*In vivo* validation of free-space fluorescence tomography using nuclear imaging,” Opt. Lett. 35(18), 3024–3026 (2010).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.35.003024Google Scholar

L. Herveet al., “Noncontact fluorescence diffuse optical tomography of heterogeneous media,” Appl. Opt. 46(22), 4896–4906 (2007).APOPAI0003-6935http://dx.doi.org/10.1364/AO.46.004896Google Scholar

G. AlexandrakisF. R. RannouA. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/50/17/021Google Scholar

V. NtziachristosR. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized born approximation,” Opt. Lett. 26(12), 893–895 (2001).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.26.000893Google Scholar

V. Ntziachristoset al., “Diffuse optical tomography of highly heterogeneous media,” IEEE Trans. Med. Imag. 20(6), 470–478 (2001).ITMID40278-0062http://dx.doi.org/10.1109/42.929613Google Scholar

R. WeisslederV. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9(1), 123–128 (2003).1078-8956http://dx.doi.org/10.1038/nm0103-123Google Scholar

A. SoubretJ. RipollV. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized born ratio,” IEEE Trans. Med. Imag. 24(10), 1377–1386 (2005).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2005.857213Google Scholar

D. Hydeet al., “Performance dependence of hybrid x-ray computed tomography/fluorescence molecular tomography on the optical forward problem,” J. Opt. Soc. Am. A 26(4), 919–923 (2009).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.26.000919Google Scholar

W.-F. CheongS. A. PrahlA. J. Welch, “A review of the optical properties of biological tissues,” J. Quant. Elect. 26(12), 2166–2185 (1990).IEJQA70018-9197http://dx.doi.org/10.1109/3.64354Google Scholar

A. B. Milsteinet al., “Fluorescence optical diffusion tomography,” Appl. Opt. 42(16), 3081–3094, (2003).APOPAI0003-6935http://dx.doi.org/10.1364/AO.42.003081Google Scholar

S. R. Arridge, “Photon-measurement density functions. part i: analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995).APOPAI0003-6935http://dx.doi.org/10.1364/AO.34.007395Google Scholar

S. R. ArridgeM. Schweiger, “Photon-measurement density functions. Part 2: finite-element-method calculations,” Appl. Opt. 34(34), 8026–8037 (1995). APOPAI0003-6935http://dx.doi.org/10.1364/AO.34.008026Google Scholar

M. Schweiger, Application of the finite element method in infrared image reconstruction of scattering media, Ph.D. Thesis, University College London (1994).Google Scholar

M. Schweigeret al., “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).MPHYA60094-2405http://dx.doi.org/10.1118/1.597634Google Scholar

R. B. Schulzet al., “Hybrid system for simultaneous fluorescence and x-ray computed tomography,” IEEE Trans. Med. Imag. 29(2), 465–473 (2010).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2009.2035310Google Scholar

X. Inteset al., “Projection access order in algebraic reconstruction technique for diffuse optical tomography,” Phys. Med. Biol. 47(1), N1–N10 (2002).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/47/1/401Google Scholar

T. StrohmerR. Vershynin, “A randomized kaczmarz algorithm with exponential convergence,” J. Fourier Anal. Appl. 15(2), 262–278 (2009).1069-5869http://dx.doi.org/10.1007/s00041-008-9030-4Google Scholar

B. Dogdaset al., “Digimouse: a 3-D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/52/3/003Google Scholar