In this paper we restrict our attention to problems in which the transport physics [

2] associated with the propagation of light at wavelength

*λ* through tissue can be adequately approximated using a diffusion model [

23,

24] of the form

where Φ(

**r**,

*λ*) is the photon fluence rate at position

**r** due to light of wavelength

*λ* injected into the medium,

*v* is the electromagnetic propagation velocity in the medium,

is the spatially varying absorption coefficient, and

*S*(

**r**,

*λ*) is the photon source with units of optical energy per unit time per unit volume. For the work in this paper the sources are considered to be delta sources in space and can be written as

*S*(

**r**,

*λ*) =

*S*_{0}(

*λ*)

*δ*(

**r** −

**r**_{s}) with

*S*_{0}(

*λ*) the source power at wavelength

*λ*. For spatially varying scattering we assume that the diffusion coefficient

*D*^{0}(

*r*,

*λ*) follows Mie scattering theory where a scattering prefactor Ψ depends on the size and density of scatterers while a scattering exponent

*b* depends on the size of the scatterers. Using this, we write the perturbation in the diffusion coefficient as

The arbitrarily chosen reference wavelength

*λ*_{0} is introduced to achieve a form of the Mie model where Ψ has the units of mm

^{−1} and Ψ′ has units of mm. In this paper, for simplicity we consider the case where the background medium is infinite and homogeneous. Generalization to the more practical case where there are boundaries is straightforward in theory though somewhat more complex in terms of implementation [

25]. As the primary objective of this work is the demonstration of a new approach to inversion, we prefer to consider the simpler physical problem holding off on the more realistic implementation for the future. Additionally, we only consider continuous wave experiments since this is what our instrument measures, where

Eq. (1) is usually written with a

*jω*/

*D*(

*λ*) term on the right hand side, in our case we consider

*ω* = 0 giving us the form shown in

Eq. (1).

We employ the Born approximation by decomposing

and

*D*^{0}(

**r**,

*λ*), as the sum of a constant background absorption,

*μ*_{a}(

*λ*), and a spatially varying perturbation Δ

*μ*_{a}(

**r**,

*λ*) as well as constant background diffusion

*D*(

*λ*) and a diffusion perturbation Δ

*D*(

**r**,

*λ*). To obtain a linear relationship between the measurements and the chromophore concentrations, we subtract

Eq. (1) from the perturbed version of the diffusion equation which gives

where

and Δ

*k*^{2}(

**r**,

*λ*) = (

*v*/

*D*(

*λ*))Δ

*μ*_{a}(

**r**,

*λ*). Assuming the availability of a Green’s function,

*G*(

**r**,

**r**′,

*λ*) for the solution of

Eq. (3) as is the case for an unbounded medium as well as range of bounded problems [

26], we can change

Eq. (3) into an integral Eq. [

27]

where

**r**_{d} is the location of the detector and (with a small abuse of notation) Φ

_{i}(

**r**,

**r**_{s},

*λ*) is used here to denote the incident field at position

**r** and wavelength

*λ* due to a delta-source located at

**r**_{s}. It should also be noted that we obtain this equation under the assumption that the total fluence rate, Φ, can be approximated as the incident field, Φ

_{i}, since Φ

_{i} Φ

_{s} [

2].

Equation (4) provides a linear relationship between the scattered fluence rate and the absorption perturbation. To relate the scattered fluence rate to concentrations of chromophores, Δ

*μ*_{a} is decomposed as follows [

28]

where

*N*_{s} is the number of absorbing species for the problem under investigation,

*ε*_{k}(

*λ*) is the extinction coefficient for the

*k*^{th} species at wavelength

*λ*, and

*c*_{k}(

**r**) is the concentration of species

*k* at location

**r**. To obtain the fully discrete form of the Born model used in Section 4, we expand each

*c*_{k}(

**r**)

where

*c*_{k,j} is the value of the concentration for species

*k* in

*V*_{j}, the

*j*^{th} “voxel”. The

(

**r**) function is an indicator function where

After using

Eqs. (5) and

(6) in

Eq. (4) and performing some algebra we obtain

We approximate

Eq. (4) as the value at the center of each pixel multiplied by the area or volume of each pixel or voxel, so in

Eq. (8) we use

*a* as the area of a pixel. This setup is illustrated in .

Setting up the linear algebraic structure associated with

Eq. (8) we define

**c**_{k}

^{Nv} as the vector obtained by lexicographically ordering the unknown concentrations associated with the

*k*^{th} chromophore and

**c**_{k+1} = ΔΨ′ and Φ

_{s}(

*λ*) is the vector of observed scattered fluence rate associated with all source-detector pairs collecting data at wavelength

*λ*. Now, with

*N*_{λ} the number of wavelengths used in a given experiment,

Eq. (8) is written in matrix-vector notation as

It should be noted in

Eq. (9) that the matrix has elements which are also the matrices

and

. The (

*m*,

*j*)

^{th} element of the

matrix is given by the product (

*v*/

*D*(

*λ*_{l}))

*G*(

**r**_{m},

**r**_{j},

*λ*_{l})Φ

_{i}(

**r**_{j},

**r**_{m},

*λ*_{l}), where

*m* represents the

*m*^{th} source-detector pair and as before

*j* represents the

*j*^{th} voxel. For

the matrix elements are given by (

*v*/

*D*(

*λ*_{l}))

*G*(

**r**_{m},

**r**_{j},

*λ*) ·

Φ

_{i}(

**r**_{j},

**r**_{m},

*λ*)Δ

*D*_{j}(

*λ*).

To evaluate the forward model for realistically sized problems, we compute the *N*_{λ} matrices **K**^{a} and **K**^{d} and store it along with the *N*_{λ} × *N*_{c} extinction coefficients. This reduces the amount of memory required for the reconstruction.