A Parallel Reconstruction Scheme in Fluorescence Tomography Based on Contrast of Independent Inversed Absorption Properties

Based on an independent forward model in fluorescent tomography, a parallel reconstructed scheme for inhomogeneous mediums with unknown absorption property is proposed in this paper. The method considers the two diffusion equations as separately describing the propagation of excited light in tissues with and without fluorescent probes inside. Then the concentration of fluorophores is obtained directly through the difference between two estimations of absorption coefficient which can be parallel inversed. In this way, the multiparameter estimation problem in fluorescent tomography is transformed into two independent single-coefficient determined schemes of diffusion optical tomography (DOT). Any algorithms proved to be efficient and effective in DOT can be directly applied here. In this study the absorption property is estimated from the independent diffusion equations by a gradient-based optimization method with finite element method (FEM) solving the forward model. Simulation results of three representative occasions show that the reconstructed method can well estimate fluorescent property and tissue absorption distribution.


INTRODUCTION
The mighty advance in biocompatible, specific fluorescent probes, combining with progress in photonic technology and methodology of modeling photon propagation, greatly promoted the development of fluorescence tomography in recent years. By transporting the ways gene expression and protein function visualized in vivo, fluorescence tomography becomes a promising means for molecularly based noninvasive imaging of biological tissues [1,2].
Always fluorescent beacons emitting in the near infrared (NIR) bandwidth are preferred, since in this spectral window hemoglobin and water absorb minimally so as to allow photons to penetrate for several centimeters in tissues. In continuous wave (CW) mode, the following coupled diffusion equations are always employed to describe the propagation of both excitation and fluorescent emission light in diffusive medium [3][4][5]: ÖD x (r)ÖΦ x (r) μ ax (r) + μ a f (r) Φ x (r) = Θ s δ r r sk , where Φ x,m is the photon density for excitation (subscript x) or fluorescent light (subscript m), D x,m (r) is the diffusion coefficient, and μ ax,m (r) is the absorption coefficient of the tissue. The absorption of excite light due to fluorophores is described as μ a f (r) and the fluorescent yield ημ a f (r) is the required fluorescence parameter. Based on this model, several simulations with 3-dimensional and 2-dimensional geometries have been done using inversion approaches such as Born approximation employing Green's function [6,7], Newton's or Newton-type optimization methods [3,4] and disturbed method [8]. However, these studies only deal with reconstruction of fluorescence properties in homogeneous mediums or under the assumption of known background optical properties. For inhomogeneous mediums with unknown absorption property as in many practical cases, only A. B. Milstein et al. reconstructed fluorescent yield and lifetime as well as optical properties in frequency mode using nonlinear Bayesian method [5]. However, in their work the fluorescent properties could only be estimated from (1b) after the optical coefficients are determined. As the two inversion parts must be executed serially, the total reconstruction process is very time consuming.
2 International Journal of Biomedical Imaging A recently developed model makes the two diffusion equations independent and speeds up the forward model execution [9]. Based on this model, we proposed a new parallel reconstructed scheme for inhomogeneous mediums with unknown absorption property. In this study the absorption property is separately inversed from the independent diffusion equations by a gradient-based optimization method with finite element method (FEM) solving the forward model. Simulation results of three different occasions show that the reconstructed method can well estimate fluorophore concentration and tissue absorption distribution.

Independent forward model
It can be easily seen that the two equations in (1) are originally dependent. For example, in (1a) the fluorescent probe acts as an absorptive heterogeneity to excited light with absorption coefficient μ a f (r), which equals ε x c(r). Here c(r) is the spatially varying fluorophore concentration and ε x is the molar extinction coefficient and is always assumed constant. In the right items of (1b) the strength of fluorescent source is influenced by the excited photon flux at that position, which mainly depends on the incident excited light and tissue's optical property at the excitation wavelength. The fluorophores absorb the excited light and then emit fluorescent light with the probability of η, the quantum yield, which depends on the type of fluorophore and the chemical environment [5].
In CW mode the diffusive coefficient is always considered spatially invariant and known in advance as D x,m . For fluorescent beacons whose excited and emission light are both in the near infrared spectral window, the Stokes shift is relatively small and the optical properties of the two wavelengths are similar, such as in human breast tissue [10,11]. Thus the optical coefficients for both wavelengths λ x and λ m can be considered equal, such as in [5,6], that is, D x = D m = D and μ ax = μ am = μ a . Based on above assumptions, R. B. Schulz et al. [9] proposed an independent forward model. Firstly they rewrite (1) as Then (2a) is converted to the following form: With respect to Φ t = (1/η)Φ m + Φ x , (4a) and (4b) are independent and can be solved simultaneously. A parallel FEM system was proposed in [9] which cuts down by one half the time of original forward problem solving.
Fluorophore concentration Parallel execution Figure 1: The process of parallel inversion scheme. The two "F" at the left are short for "fluorophore." The reconstruction of fluorescence tomography in heterogeneous medium is converted to two parallel absorption coefficient inversed problems in DOT.

Parallel inversion scheme
It can be easily found that the independent equations (4a) and (4b) have totally the same form as the single diffusion equation in diffuse optical tomography (DOT). The source items on the right side of equation are both the same δ functions and all the unknown properties are moved into the coefficient positions of the partial differential equations. Thus, for each equation, the estimation of unknown properties becomes a nonlinear reconstruction problem right as in DOT. Equation (4) also has its physical significance, although educed through pure mathematic transform. Firstly (4a) describes the propagation of excitation light with the absorption of both tissue itself and the inside fluorophores. Then let us focus on (4b). The quantum yield, η, is always defined as a ratio between emitting fluorescence photon numbers and the number of excitation photons absorbed by fluorophore. Therefore, the item (1/η)Φ m (r) can be treated as compensating the excitation photon density absorbed by fluorescent beacons and Φ t shows the total excitation photon flux together with the fluorophore absorbed part. Thus, (4b) describes the transportation of excitation light in tissue with assumed no fluorochrome inside, providing a referenced absorption situation of only the tissue originating absorption μ a (r) without fluorescent beacon administration.
Based on above analysis, (4a) and (4b) provide two different absorption situations of excitation wavelength with contrast caused by fluorescent beacon's absorption ε x c(r). Knowing that (4a) and (4b) have entirely the same form as absorption parameter inversed problems in DOT, we can parallel reconstruct μ a (r) + εc(r) and μ a (r) from the boundary measurement of excited light and fluorescent light according to the two independent diffusion equations. Then their contrast caused by fluorophore ε x c(r), which is proportional to the concentration of fluorescent probes, can be directly obtained by the difference of two absorption images. Figure 1 shows the process of the parallel reconstruction scheme.
So the parallel scheme converts the reconstruction of fluorescent property in heterogeneous medium to two independent absorption coefficient inversion problems in DOT. Both the unknown tissue absorption property and fluorophore concentration are determined. In this way, the multiparameter estimation problem in FODT is simplified to two single coefficient determined schemes of DOT. Any algorithms Xiaolei Song et al. proved to be efficient and effective in DOT can be directly applied here. Additionally, the total reconstructed time is the same as single DOT problem for which the two absorption property inversions can be performed at the same time.

Gradient-based optimization inversed method
To demonstrate the parallel reconstruction strategy, in this work a gradient-based optimization inversion method is used for the absorption coefficient inversion with FEM solving the forward model [12,13]. Considering an experimental setting that includes S point excitation light sources located at ξ j ¾ ∂ Ω ( j = 1, . . . , S), and M j measurement positions σ j,i ¾ ∂ Ω (i = 1, . . ., M j ) for each source j, the following objective function can be defined: where Γ j,i represents the photon intensity measured at position σ j,i with the incident excitation source located at ξ j . The subscript c denotes the values calculated by the forward simulation and me represents the experimental (or in this case computer-simulated) values. By minimizing the objective function, we can handle an iteration computing and reconstruct the distribution of absorption coefficient.
To solve the optimization problem using a conjugate gradient (CG) method, the gradient of the objective function needs to be calculated as follows: Therefore, the gradient vector z can be presented as where J is an M TOT ¢N TOT Jacobian matrix, M TOT = S j=1 M j is the total measurement number at the boundary, and N TOT is the number of the coefficients to be reconstructed. Here b is the residual error between the boundary measurements and computation values. Then J can be calculated by an adjoint source scheme based on the establishment of PMDF (photon measurement density function, as defined in [14]).
With the gradient calculated, the next step is to conduct one-dimensional search in order to find a best step length on this gradient direction. Then, we refresh the absorption coefficient and recalculate the gradient to form iteration computing until the error reaches the supposed value.

SIMULATION AND RESULTS
In this section, a group of simulation tests is done to demonstrate our parallel inversion scheme. For each reconstruction, three kinds of figures are produced, they are (1) absorption distribution due to both medium heterogeneity and extrinsic fluorescent probes, that is, μ a (r) + εc(r); (2) absorption coefficient distribution μ a (r) of the heterogeneous medium without fluorophore; (3) εc(r), which is just the minus result of the former two absorption coefficient distributions, representing the fluorochrome concentration as ε is always assumed to be constant.
In our simulation, the imaging area is a circle of 32 mm diameter with the Robin boundary condition Φ(ξ) + 2D ¡ AnÖΦ(ξ) = 0, ξ ¾ ∂Ω, applied [15]. For both the excitation and emission wavelengths, the diffusion parameter of the medium is assumed constant and known as  ( 6,4) and is 0.025 mm 1 in the rest. εc(r) is 0.075 mm 1 in the right 8 mm-diameter circle centred at (4,6) and is zero in the rest.
with the reduced scattering coefficient μ ¼ s = 2.0 mm 1 and the homogenous background absorption coefficient μ a = 0.025 mm 1 . The 32 excited sources and 32 detectors are positioned alternately around the circle with equal intervals between each other. As the sources are collimated beams, they are represented by isotropic point sources located at one diffusion length of 1/μ ¼ s below the border [15]. The forward model was solved by the FEMLAB2.1 (COMSOL Inc., Burlington, MA 01803, USA, a FEM solution in MATLAB (The Math Works, Inc.)), and all the codes are programmed in MATLAB6.0. The finite-element mesh for reconstruction consists of 544 nodes and 966 triangle elements (Figure 2(a)). The final estimated images reported in the paper are all of 25 CG iterations as the objective function descends very slowly when the iterative time is more than 25 (Figure 2(b)).
The simulated boundary measurement data Φ x (ξ) and Φ m (ξ) were generated separately from the forward model of (4a) and (4b) with exact distributions of μ a (r) and εc(r). According to different positions of the fluorochrome with respect to original absorption heterogeneity of tissue, we reconstructed the following three kinds of occasions in noisefree environment. (Here we abbreviate the original absorption heterogeneity as HA and absorption due to fluorescent probes as FA), that is, 1.HA and FA are apart from each other; 2.HA and FA intersect; 3.FA is within HA. We will present these results one by one.  (6,2) and is 0.025 mm 1 in the rest. εc(r) is 0.05 mm 1 in the upper 8 mm-diameter circle centred at (6,2) and is zero in the rest. Thus, the total absorption in the overlapped part is 0.10 mm 1 in (a). reconstructed qualities of the two dependent inverse problems, the minus result (c) will have values less than zero. However it is known that the concentration should not be negative in practice. So that these negative values are constrained to zero in our simulation. The modified results with nonnegative restriction in (d) show less background noise and better reconstructed qualities compared with the unconstrained ones of (c) in each figure.

DISCUSSION AND CONCLUSIONS
In last section, through the parallel inversion scheme, we reconstructed the three different relative distributions of tissue absorption heterogeneity and fluorophore in 2-dimensional geometry, combined with a nonnegative restriction. Distri-butions of fluorophore concentration as well as the absorption properties of tissue is obtained, by treating fluorescent tomography in inhomogeneous tissues as two independent nonlinear estimations in DOT. Although reconstruction of fluorophore concentration is a linear problem when the optical properties are known in advance, the parallel method has its superiorities. Firstly the two nonlinear estimations can be executed synchronously, so the total reconstruction time keeps the same as single absorption inversion in DOT. Then, any reconstruction algorithms which are already proved to be effective and efficient in DOT can be directly applied here. In this paper a basic CG method is employed for the reconstruction of absorption coefficient with noise-free data. However, for noisy data, the spatial resolution and the ability to solve deep located heterogeneities are very limited with only the  Figure 5: Images inversed when FA and HA intersect. (a) and (b) are, respectively, the independent reconstructed images of μ a (r)+εc(r) and μ a (r). (a) minus (b) educed (c), describing the concentration of fluorochome. (d) constrained the negative values in (c) to zero. The black circles represent the original exact locations and shapes of HA and FA. The original μ a (r) is 0.05 mm 1 in the bigger 12 mm-diameter circle centred at (4, 0) and is 0.025 mm 1 in the rest. εc(r) is 0.05 mm 1 in the inside 5 mm-diameter circle centred at (6, 0) and is zero in the rest. basic CG method [13]. Other strategies can be added to improve the reconstructed image quality such as the spatial located weighted method [13] and adaptive finite element technology [16]. Additionally, it does not require measurements of excitation light before fluorochrome administration for estimating the absorption property of tissue, as Φ t in (4b) could provide an absolute photon-field measurement.
In the independent model and the parellel reconstruciton method, a constant quantum yield is a basic assumption of quantitative fluorescence imaging as in [5,17]. In fact, it might be questionable in biological environments. However, newly developed emitters like quantum dots which are embedded in an isolator are independent of the environment [18]. Therefore, it offers a practical and relative simple method to determine fluorescent concentration in heterogeneous medium and tissues, as Φ t can be easily combined from boundary measurements Φ x , Φ m using the optics filter.