Time-gated perturbation Monte Carlo for whole body functional imaging in small animals

Jin Chen; Xavier Intes

doi:10.1364/OE.17.019566

1. Introduction

Preclinical or small-animal in-vivo imaging has experienced an exponential phase of growth in the last decade [1, 2]. In vivo pre-clinical imaging enables researchers to identify biologic processes, monitor the efficacy of compounds and measure the effects of disease progression over time in intact host environments. It is becoming an essential translational tool between in vitro research and clinical application.

The blooming use of optical techniques for in vitro and in vivo research has generated great interest in developing new imaging platforms, with clear potential for translation into the clinic [3]. As a non-linear process dominated by scattering, light propagation in tissue is highly dependent on the optical properties of the sample imaged [4]. Thus, optical imaging technique requires accurate quantification of the absolute optical properties of the probed specimen, which is essential for rigorous photon propagation modeling in tissue and, hence, important for improving optical reconstruction fidelity [5, 6].

Diffuse Optical Tomography (DOT) is a well established technique enabling thick-tissue three-dimensional imaging that offers the unique benefit of direct quantitative recovery of functional tissue properties [7, 8]. The technique relies on sequentially exciting multiple points on the tissue boundary and on collecting diffuse light exiting a few centimeters away from the excitation points with multiple detectors. Combined with sophisticated illumination and detection schemes, physical or mathematical models of photon propagation in tissue are employed to quantitatively interrogate the volume to be imaged. The technique has been explored broadly in clinical settings and has been applied successfully to small animal fluorescence imaging [9].

However, pre-clinical implementation of DOT faces several unique challenges that need to be addressed in order to produce robust and truly quantitative imaging platforms for wide-spread acceptance of this new modality. First, data type selection plays an important role in the DOT technique. Three optical technology domains exist in the diffuse optical imaging field. Continuous wave (CW) is the most wide-spread technique due to its simplicity to implement and low cost components; however it provides limited information content leading to non-unique solution in the inverse problem [10]. Frequency domain (FD) measures modulation amplitude attenuation and phase shift due to light propagation which allows one to separate the contributions of absorption and scattering [11, 12]. However, producing instruments with dense source- detectors arrays working at multiple frequencies in the range of a few MHz to 1 GHz is extremely challenging [13]. Time-domain (TD) methods illuminate tissue with ultrafast photon pulses and resolve the arrival of the photons as a function of time at different locations around the tissue boundary. The temporal distribution of photons detected after propagation is known as the temporal point spread function (TPSF). TD technology offers a rich data set that is unmatched by CW and FD methods for quantitative accuracy, minimal cross talk between chromophores, and improved resolution [14–17]. Furthermore, TD systems are more sensitive than their CW counterparts [18].

The second challenge facing pre-clinical DOT is the structural complexity of the volume to be imaged. Whole body small animal imaging is characterized by complex boundaries and complex spatial distribution of small organs that exhibit a wide range of optical properties. Employing forward models based on the diffusion approximation in such highly heterogeneous medium, where low scattering and high absorption regions exist, leads to erroneous data interpretation [19, 20]. The perturbation Monte Carlo (pMC) method is proposed as an alternative to the diffusion approximation based tomographic formulation [21]. The Monte Carlo technique is easy to implement, flexible, proven to be accurate over a large span of optical properties (or conversely over a large spectral range), accurate for small interrogated volumes and appropriate for complex boundaries. Since it is a statistical method, a large number of packets of photons need to be simulated to obtain simulations with adequate signal to noise ratio. This translates into large computing power requirements to achieve acceptable simulations within hours. This drawback has confined the Monte Carlo method to simulate light propagation as a validating tool for other forward models or to be used to estimate bulk optical properties of tissues. However, thanks to recent advances in parallel computing, optical tomography of small animals based on a computationally efficient parallel Monte Carlo forward model has lately been achieved in the time domain [22].

The third challenge lies in the ill-posed nature of the inverse problem leading to poor spatial resolution and reconstructions that are highly sensitive to noise. This ill-posedness can be mitigated by a priori knowledge of the spectral and structural properties of the imaged sample [12]. Many strategies capitalizing on anatomical priors from a high resolution imaging modality such as MRI or X-ray have been devised for optical tomography to improve the recovery of optical features and contrast [23], with methods that combine both anatomical and functional priors to guide the optical reconstruction offering superior quantitative performances for functional imaging applications [24].

Towards high resolution whole body quantitative functional imaging of small animals, we explore the combination of new optical data types, rigorous photon propagation model and anatomy guided inverse problem formulation. In this in silico investigative work, a multispectral time-gated perturbation Monte Carlo method and a hierarchical Bayesian reconstruction are applied on an anatomically accurate mouse model. This paper is structured as follows: Section 2 briefly describes the forward Monte Carlo model, the direct functional imaging method and the implementation of the Bayesian method. Section 3 describes the mouse phantom and outline of the simulation configurations. Section 4 reports the results on the computational efficiency, the sensitivity map computations, and the simulated image reconstruction and accuracy gains of the proposed method on the mouse model. Section 5 summarizes the innovations and describes the ongoing work to implement this new approach.

2. Methodology

2.1 Time-resolved perturbation Monte Carlo

The perturbation Monte Carlo (pMC) method provides an efficient framework to perform optical tomography based on the MC forward model. In this work we extend the pMC approach to construct time-resolved Jacobian for time-gated diffuse optical tomography. A comprehensive description of the photon propagation rules in MC modeling can be found in previous studies [25–27].

In the MC simulations, with a set of initial optical properties, the unperturbed measurement w of one photon packet at a detector can be obtained. If the optical properties change in a certain region (perturbed region) along the photon path, w should be modified to a perturbed measurement $\hat{w}$ . A prediction for $\hat{w}$ of this perturbed photon packet at the detector can be estimated by:

\hat{w} = w {(\frac{{\hat{μ}}_{s} / {\hat{μ}}_{t}}{μ_{s} / μ_{t}})}^{p} {(\frac{{\hat{μ}}_{t}}{μ_{t}})}^{p} \exp (- ({\hat{μ}}_{t} - μ_{t}) L) .

where

{\hat{μ}}_{s} = {\hat{μ}}_{s} + δ μ_{s}

,

{\hat{μ}}_{a} = μ_{a} + δ μ_{a}

,

μ_{t} = μ_{a} + μ_{s}

and

{\hat{μ}}_{t} = {\hat{μ}}_{a} + {\hat{μ}}_{s}

, with

μ_{a}

and

μ_{s}

being the initial values of the absorption and scattering coefficient, P is the number of scattering events and L is the path length in the perturbed region with absorption change

δ μ_{a}

and scattering change

δ μ_{s}

. This equation provides an accurate means to calculate the weight for every single photon exiting the sample and efficiently compute the absorption and scattering Jacobians [21, 28] to perform DOT.

It is straightforward to extend the pMC method to time resolved studies. Assuming that the perturbations are time invariant over the propagation time (~ns), the time dependent photon weight becomes:

\hat{w} (t) = w (t) {(\frac{{\hat{μ}}_{s} / {\hat{μ}}_{t}}{{\hat{μ}}_{s} / {\hat{μ}}_{t}})}^{p} {(\frac{{\hat{μ}}_{t}}{μ_{t}})}^{p} \exp (- ({\hat{μ}}_{t} - {\hat{μ}}_{t}) L) .

where

\hat{w} (t)

and

w (t)

are the photon weights detected in a time-bin. When time-gated pMC is applied to voxelized three dimensional imaging, each voxel is considered to be an independent perturbed region. Therefore storing all photon paths as in the classic pMC becomes inextricable. To cope with the large number of unknown regions, we collect the photons exiting at each detector position for each time-gate. We store the accumulated detector readings and the weighted average number of collisions and length at each voxel. The final estimation of the accumulated weight

\hat{W} (t)

can then be expressed as:

\hat{W} (t) = W (t) {\prod_{j = 1}^{n} (\frac{{\hat{μ}}_{s} (r_{j}) / {\hat{μ}}_{t} (r_{j})}{μ_{s} (r_{j}) / μ_{t} (r_{j})})}^{p (r_{j}, t)} {(\frac{{\hat{μ}}_{t} (r_{j})}{μ_{t} (r_{j})})}^{p (r_{j}, t)} \exp (- ({\hat{μ}}_{t} (r_{j}) - μ_{t} (r_{j})) L (r_{j}, t)) .

where

W (t)

is the unperturbed detector reading,

P (r_{j}, t)

is the number of collisions and

L (r_{j}, t)

is the path length in voxel

r_{j}

at time gate t. n is the total number of the voxels. With

δ μ (r_{j})

the perturbation in optical coefficients (either

μ_{a}

or

μ_{s}

) at

r_{j}

, Eq. (3) provides a means to estimate the changes in photon weight introduced by changes in

δ μ (r_{j})

and is then used to find the terms of the Jacobian matrix

J (t)

with

J_{i j} (t) = \partial {\hat{W}}_{i} (t) / \partial δ μ (r_{j})

for this time gate. At a specific time-gate t, the problem can be then casted in the formulation:

[\begin{array}{c} Δ W_{1} (s d, t) \\ ⋮ \\ Δ W_{m} (s d, t) \end{array}] = [\begin{array}{c} J_{11} (t) & \dots & J_{1 n} (t) \\ ⋮ & ⋱ & ⋮ \\ J_{1 m} (t) & \dots & J_{m n} (t) \end{array}] \times [\begin{array}{c} Δ δ μ (r_{1}) \\ ⋮ \\ Δ δ μ (r_{n}) \end{array}] = J (t) \times [\begin{array}{c} Δ δ μ (r_{1}) \\ ⋮ \\ Δ δ μ (r_{n}) \end{array}] .

where

Δ W_{i} (s d, t) = {\hat{W}}_{i} (s d, t) - W_{i} (s d, t)

is the total photon weight difference for the

i^{t h}

selection of source-detector pair

s d

and time-gate t and

J (t)

is the Jacobian matrix. This equation can be easily expanded for the multiple gates case:

[\begin{array}{c} Δ W_{1} (s d, t_{1}) \\ ⋮ \\ Δ W_{m \times k} (s d, t_{k}) \end{array}] = [\begin{array}{c} J (t_{1}) \\ ⋮ \\ J (t_{k}) \end{array}] \times [\begin{array}{c} Δ δ μ (r_{1}) \\ ⋮ \\ Δ δ μ (r_{n}) \end{array}] = J \times [\begin{array}{c} Δ δ μ (r_{1}) \\ ⋮ \\ Δ δ μ (r_{1}) \end{array}] .

where k is the number of selected gates and

m \times k

is the total number of selections of measurements. According to the method discussed above, two sets of Monte Carlo simulations are required for in silico investigations; one is to perform a simulation with the initial guess (or background) of optical properties. It is worthy to note that the initial guess does not have to be homogeneous; any heterogeneity in the phantom model can be considered. This simulation generates the unperturbed detector readings W and the Jacobian J. The other set is to compute the perturbed measurement

\hat{W}

with changes in optical properties. A linear system can be formed after these two simulations are completed. Such an approach has the benefit to provide a set of measurements that are based on the most rigorous forward model of light propagation.

2.2 Functional imaging

Assuming that deoxy-hemoglobin (Hb), oxy-hemoglobin(HbO₂) and water (H₂O) are the only chromophores contributing to the global absorption in the spectral window considered, the physiological information or the changes in the concentration of Hb, HbO₂ and H₂O at different wavelengths for $j^{t h}$ voxel are linearly related to the absorption coefficients as follows:

δ μ_{a}^{λ_{k}} (r_{j}) = [ε_{H b}^{λ_{k}} ε_{H b O_{2}}^{λ_{k}} ε_{H b O_{2}}^{λ_{k}}] \times [\begin{array}{l} δ [H b (r_{j})] \\ δ [H b O_{2} (r_{j})] \\ δ [H_{2} O (r_{j})] \end{array}]

where

δ μ_{a}^{λ_{k}} (r_{j})

is the differential absorption coefficient at the

k^{t h}

wavelength

λ_{k}

and for the

j^{t h}

voxel,

ε_{H b}^{λ_{i}}

,

ε_{H b O_{2}}^{λ_{i}}

and

ε_{H_{2} O}^{λ_{i}}

are the extinction coefficient of the corresponding chromophores at the

k^{t h}

wavelength. Combining this relationship to Eq. (6), we obtain:

[\begin{array}{c} δ W_{1}^{λ_{k}} (s d, t) \\ ⋮ \\ δ W_{m}^{λ_{k}} (s d, t) \end{array}] = [ε_{H b}^{λ_{k}} J^{λ_{k}} ε_{H b O_{2}}^{λ_{k}} J^{λ_{k}} ε_{H_{2} O}^{λ 1} J^{λ_{k}}] \times [\begin{array}{c} δ [H b] \\ δ [H b O_{2}] \\ δ [H_{2} O] \end{array}]

where

δ [H b] = [δ [H b (r_{1})] \dots δ [H b (r_{n})]]'

,

δ [H b O_{2}] = [δ [H b O_{2} (r_{1})] \dots δ [H b O_{2} (r_{n})]]'

and

δ [H_{2} O] = [δ [H_{2} O (r_{1})] \dots δ [H_{2} O (r_{n})]]'

.

δ W_{1}^{λ_{k}} (s d, t)

and

J^{λ_{k}}

are the differential detector reading and the Jacobian at the

k^{t h}

wavelength

λ_{k}

, respectively. Multispectral priors can thus be incorporated into the tomographic inverse problem to directly retrieve functional maps [29].

2.3 Inverse problem

To overcome the ill-posedness associated with whole body imaging, we employ a Bayesian formulation to incorporate anatomical and functional priors. This formulation constrains the space of unknowns with a spatially varying probability density function derived from high-resolution anatomical maps and physiological priors that are implicitly defined to constrain the reconstructions in physiologically meaningful ranges. The full derivation of the theoretical developments for the inverse formulation and algorithm is presented in detail in references [24, 30]. A brief description of the method is provided below.

Having designed the hierarchical noise and image models, we formulate the joint distribution of the measurements, the image and the unknown hyper-parameters associated with the noise and image models. In order to estimate the hyper-parameters, we use a conjugate gradient (CG) algorithm for a hyper-parameter estimation step followed by an image update. The maximum a posteriori (MAP) estimates of the hyper-parameters at each CG iteration prior to the image update is:

{\hat{x}}_{M A P} = \arg \max_{x} {\log p (y | x) + \log p (x, σ | C)}

where y is the measurement vector,

p (y | x)

is the data likelihood, σ is the vector of variance of sub-images and

p (x, σ | C)

is the conditional probability prior on the unknown image x and variance σ given the tissue label information C. The label information C can be taken into account by considering one sub-image with the same label at a time. The optimization problem can be solved by alternatively updating these two parameters image x and the variance σ. Assuming a Gaussian distribution for x, the probability density function of

i^{t h}

sub-image

x_{i}

is:

p (x_{i} | σ_{i}) = \frac{1}{{(2 π σ_{i}^{2})}^{N_{i} / 2}} \exp (- \frac{1}{2 σ_{i}^{2}} | | x_{i} - μ_{i} | |^{2}), i = 1, 2, \dots, M

where M is the number of sub-regions,

N_{i}

is the number of voxels in the

i^{t h}

sub-image, and

μ_{i}

is the mean value of this sub-image, Similarly, assuming a Gaussian distribution for the standard deviation

σ_{i}

of the

i^{t h}

sub-image, we have:

p (σ_{i}) = \frac{1}{{(2 π γ_{i}^{2})}^{N_{i} / 2}} \exp (- \frac{1}{2 γ_{i}^{2}} | | σ_{i} - {\bar{σ}}_{i} | |^{2}), i = 1, 2 \dots, M

where

γ_{i}

is the variance and

{\bar{σ}}_{i}

is the mean value of variance

σ_{i}

.

γ_{i}

and

{\bar{σ}}_{i}

are predefined as the physiological prior allowing to constrain the reconstruction. The image and noise are accommodated by respectively applying Eq. (9) and (10) for all sub-images at each update along with the solution process.

3. Simulation configuration

3.1 Geometric model and optode settings

We tested our newly developed pMC scheme using an anatomically accurate 3D mouse atlas created from CT and cryosection slices at University of Southern California [31, 32]. The mouse model is discretized to $91 \times 35 \times 22$ voxels with size 1mm³. Simulations were performed under a transmittance geometry mimicking the configuration of our time resolved multi-spectral non-contact imaging platform [33].

3.1.1 Computational efficiency

We evaluated the computational efficiency of the time-gated forward MC versus the number of CPUs, number of gates and number of detectors used in the simulations. The motivation of this evaluation is due to (a) the low computational efficiency of MC method and (b) the temporally and spatially dense data sets required in time-resolved whole body imaging.

First, we tested the computational feasibility of the proposed method using different numbers of CPUs for one source-detector pair. For this case, the source and detector were facing each other in a transmittance geometry (2cm separation), probing an homogenous media. Second, we evaluated the effect on time cost caused by changes in the number of detector and gate simulated. The recent developments in ultra-fast gated CCD cameras allow acquiring spatially dense time-resolved data sets. Such sub-millimeter spaced arrays of sources and detectors (that is, data to the order of $10^{4}$ - $10^{6}$ measurements or more) are required for high-fidelity, small animal imaging [18]. In the pMC simulations, multiple detector readings and the corresponding Jacobian can be obtained simultaneously for all gates. We investigated the time cost associated with simulating 16, 64, 256 or 1024 evenly spaced detectors. Sources and detectors were simulated with a 1mm diameter in all simulations herein. The area covered by the different detector sets is fixed at 16mm × 16mm. Additionally, for each of these source-detector settings, time-gated Jacobians for 1, 10, 20, 30, 40, 50 gates were recorded and the time cost for each simulation was evaluated.

3.1.2 Anatomically accurate simulations

To assess the performances of time-gated functional imaging based on pMC formulation for whole-body small animal imaging, we generated a synthetic mouse model bearing five organs (heart, stomach, liver, kidneys and non-inflated lungs). The original high-resolution and derived synthetic phantom employed in this work are provided in Fig. 1 . We limited the volume to be imaged to the upper torso-mid section (including kidneys), leading to a $33 \times 41 \times 20$ voxels discritized volume to reconstruct, with 1mm³ voxels.

Fig. 1 (a) The original high resolution mouse model (b) The resized mouse model with five organs (heart {teal}, stomach {green}, liver {yellow}, kidneys {red} and lungs {purple} extracted from Fig. 1(a))

Download Full Size | PDF

Due to memory limitations associated with the significant size of the linear system to invert, we selected a subset of 100 sources and 144 detectors for all reconstructions herein. We investigated two heuristic strategies to select sources and detectors positions. First, we selected a linear distribution of sources and detectors that covered all the volume to be imaged. In this configuration, sources were evenly spaced at 4mm intervals whereas the detectors were evenly separated by 3.6mm (cf Fig. 4a ). Such homogenous distribution is the optimal configuration according to Singular Value Anlysis (SVA) [11].

Fig. 4 Anatomically guided nonlinear sampling: (a) Evenly spanned source (black dots) and detectors (red dots) configuration and reconstructed (c) blood volume and (d) oxygen saturation; (b) dense source and detector pairs around the overlapping organs and the reconstructed (e) blood volume and (f) oxygen saturation. The heart and stomach have been delineated with white lines.

Download Full Size | PDF

Second, the same number of sources-detectors were selected but with a non-linear spatial sampling distribution. Guided by the structural prior, denser source-detector pairs were employed around the area of significant overlap between the (i) lungs and the heart and (ii) the liver and the stomach,. The denser source-to-source separation was 2mm and detector-to-detector separation was 1.8mm while in other area the separation of sources or detectors was doubled (cf. Fig. 4b).

3.2 Computational settings

The computations are distributively performed on a supercomputer BlueGene (CCNI at RPI) with 16384 nodes (dual CPU at 700MHz). For all simulations herein, the forward computations were running on 1024 nodes. Based on previous work [22], we launched 10⁹ photons per source. This is the optimal number of photons considering time-cost and adequate photon statistics to reliably compute time resolved Jacobians in mouse model.

Photon profiles were computed over a 1ns time period. Time gates were simulated with a 200ps width and with 20ps gate delay, leading to overall 50 simulated time gates per detector. Moreover, simulations were performed at six wavelengths evenly spanned from 700nm to 950 nm.

The inverse calculations were executed on a personal computer (i7 3GHz, 8Gb of RAM), under matlab, using a Polack Ribiere Conjugate Gradient algorithm developed in house. Due to memory limitations, only three gates were employed for the reconstructions. From the time gate profiles (TPSFs), the time-gate with the maximum photon count is selected. The gates having photon count of half the maximum value before and after the maximum gate are selected as well (cf. Fig. 3 ).

Fig. 3 Typical Jacobians for three different time gates in a mouse model. First column: TPSF simulated and selected time gates (first half maximum, maximum and last half maximum of TPSFs). Second column: Sagittal plane of the synthetic phantom and normalized Jacobians corresponding to each selected time gate. Third column: Transverse plane of the murine model with associated normalized Jacobian. The Jacobians are plotted in log scale to provide better visualization.

Download Full Size | PDF

3.3 Functional and optical parameters

The in vivo functional properties of the murine organs are still unknown to date. Average optical properties of small animals have been reported in transmittance but in vivo organ specific values have not been reported yet due to the difficulty to resolve locally such parameters. To model realistic optical properties, we simulated optical properties based on functional parameters derived from the litterature [34, 35]. Note that these values are obtained mainly from ex-vivo measurements and we expect significant deviations from in vivo values.

Although the initial guess of optical properties can be heterogeneous in the unperturbed simulations, we considered homogeneous initial optical properties in the unperturbed simulations. That is, all organs are assigned identical optical parameters as the background. This approach is consistent with the normalized Born or Rytov approaches commonly employed in the field.

For simplicity, in this investigative study, we considered the scattering parameter spatially homogenous in both unperturbed and perturbed simulations, but spectrally following a power law $μ_{s}' = A λ^{- S P}$ ( $A = 10600$ and $S P = 1.43$ ) [36]. For the volume considered herein and with the transmittance geometry, no diffeerences were found between simulations based on the anisotropy factor $g = 0.9$ or $g = 0$ . Thus g was kept constant at 0 over the full body to reduce the computational burden in the Monte Carlo simulations. The refractive index of the tissue was set at 1.37. The functional parameters that were simulated for each organs are provided in Table 1 . These functional parameters were used to calculate the corresponding optical absorption values of each organ at the six wavelenghts simulated. The correspoding values are tabulated in Table 2 .

Table 1. Functional parameters from the literature.

View Table | View all tables in this article

Table 2. Absorption coefficients for the spectral set investigation.

View Table | View all tables in this article

4 Results

4.1 Computational feasibility

Fig. 2a depicts the time cost using different nodes. On a single CPU, a one source-one detector time-resolved Monte Carlo simulation using $10^{9}$ photons takes 62 hours for completion while it can be completed under 5 minutes if 1024 nodes are used on BlueGene. It clearly shows the dramatic increase of computational efficiency thanks to parallel computing techniques.

Fig. 2 (a) Time cost with increasing number of nodes, from 10⁵ to 10¹⁰ with homogeneous optical properties. (b) Computational time when simulating different gates and detectors.

Download Full Size | PDF

Fig. 2b shows the time cost of simulations using different numbers of detectors and gates. It is worthwhile to note that the time cost increases non-linearly as the gate number increases, which indicates that the proposed approach allows one to simulate large number of gates without significant time cost penalty. Moreover, although parallel data transfer occurs in parallel computing, the time shift of using different number detectors in the forward model is within 10%, which allows to simulate dense detectors simultaneously. These two particulars higlight the potential of time-gated pMC modeling for processing dense spatial and temporal data sets such as those obtained with a time-gated CCD camera.

4.2 Time-gated reconstruction

As an example, we provide the time-resolved Jacobian for the center source-detector pair in Fig. 3. As expected, the Jacobian shape expands as time increases. Photons collected shortly after the pulse experience significantly fewer scattering events and probe a small volume, whereas late photons sample larger volumes. This difference in probing volumes provides unique information for reconstruction and is at the origin of the superior performances of time-gated data types compared to CW data types for DOT [37].

Functional direct reconstructions were performed on the synthetic mouse model bearing six organs for a variety of scenarios using six wavelengths. Due to the high dynamical range encountered in the different projections computed, we employed a SNR cut-off filter on the measurements selected to compute the inverse problem. Based on heuristic tests, all gates below 500 counts were not considered.

We investigated the effects of data type (CW or TG), source-detector sampling and soft prior constraints on the accuracy of DOT to retrieve the quantitative biodistribution of the main functional parameters simulated. The paraneters that are reported in this work are the blood volume (BV), the relative oxygen saturation (StO₂) and the water content (H₂O). The relative errors for the $i^{t h}$ sub-image associated with these parameters are defined by:

{error}_{X_{_{i}}} = ({\bar{\hat{X}}}_{i} - X_{i}) / X_{i}

where

{\bar{\hat{X}}}_{i}

is the mean of the reconstructed image over the

i^{t h}

sub-region and

X_{i}

the expected value in this sub-region (X being BV or H₂O respectively). As the StO₂ is already a percentage, we provide here the absolute error:

{error}_{X_{i}} = {\bar{\hat{X}}}_{i} - X_{i}

for X being StO₂. The standard deviation for the image of the

i^{t h}

region are calculated as

σ_{i} = \sqrt{\frac{1}{N_{i}} | | {\hat{X}}_{i} - {\bar{\hat{X}}}_{i} | |^{2}}

where

N_{i}

is the number of voxels and

{\hat{X}}_{i}

is the reconstructed image in the

i^{t h}

sub-image. These error and standard deviation functions were computed over the reference volume of each organ derived from the mouse atlas. The values of these error functions for all cases investigated are provided in Table 3 .

Table 3. Errors and standard deviations of the reconstructed maps: a crossed cell in n-Lin or Baye. indicates reconstructions using non-linear source detector sampling and Bayesian formulation in the inverse problem conversely to linear sampling and non-constrained reconstructions if not marked. All the values herein are provided in %.

View Table | View all tables in this article

An example of functional reconstructions for the case of linear and non-liner source-detector sampling with a priori constraints is provided in Fig. 4.

In all cases investigated herein, the reconstruction using time-gated data with anatomical a priori presents overall higher accuracy. Especially, combining TG and spatial constraints provides better estimates than unconstrained CW or TG reconstructions (even though spectral a priori was used). Without spatial priors, TG data outperforms CW data, leading to more accurate mean estimates, but suffers from large quantitative spatial deviations within the volume of the organs. This is in good agreement with previous results reported in reference [20].

A significant improvement in the quantitative estimation of the blood volume and the saturation was observed when non-linear dense source patterns were employed. As seen in Fig. 4, even source distributions over the whole body lead to under-sampling of important small heterogeneous volume. In this case, the non-linear pattern led to accurate parameters estimation (within 6% for each organ) whereas linearly distributed patterns did not resolve the heart or the stomach despite the soft spatial constraint.

5. Discussion

In this work, simulations were performed to asses the potential of time-resolved functional DOT for whole body small animal imaging. To tackle the challenges of small animal imaging, we explored in silico the feasiblity to combine Monte Carlo based reconstructions, time-gated data type and soft a priori constraints to quantitatively image the endogenous functional properties of murine model organs.

First, we demonstrated that the perturbative Monte Carlo formulation is suitable and feasible for functional DOT in preclinical scenarios, where the specimens present a small optically complex model with multiple close organs of significantly low-scattering (bladder and lungs) or highly absorbing tissues. The method reconstructs wide range of coefficients in tissue accurately and provides satisfactory results synthetically. While only functional parameters related to absorption coefficients are reconstructed, perturbations of structural parameters related to scattering coefficient can be addressed [37]. We also argue that the time cost for the reconstruction is acceptable for in vivo imaging applications. The low time efficiency limitation in Monte Carlo simulations is resolved by the use of massively parallel computing. We plan to further decrease computational cost by investigating other method or hardware to reduce the timecost of the forward model computations [38].

Second, we investigate a new data type for functional imaging: time gates. Such data type provides a unique set of information to perform quantitative DOT. Thanks to the different volumes probed by different time gates, Time-Gated DOT offers superior quantification and resolution compared to CW reconstructions. Time-resolved technology is a promising modality with strength relying on its ability to estimate absolute characterization of tissue optical properties as well as the structural information more precisely. Especially, time-gated technology permits one to acquire both spatially and temporally dense data sets for high-fidelity three-dimensional imaging. However, careful instrumental calibration and normalization should be applied when experimental data are employed to reduce the temporal errors. We are currently investigating new inverse formulations to reduce this experimental drawback. Also, all reconstructions herein have been performed based on a subjective gate selection. We plan to futher explore the optimal selection of time gates, especially in the context of functional and structural parameters cross-talk.

Last, incorporating a-priori anatomical constraints provided stable reconstructions and improved quantification. The 3D optical parameters estimation suffers from low resolution, partial volume effect and inter-parameter cross-talk due to the ill-posed nature of DOT. Casting spectral and structural priors into a Bayesian framework yields more accurate and stable calculations to provide enhanced in-vivo functional estimates. Without anatomical guidance, DOT is unable to recover accurately the functional properties of all the organs probed despite spectral priors and dense spatial and temporal data sets. In addition, the present article provides the first insights on the importance of anatomically guided probe settings.

The reconstruction of high absorbing and overlapping tissues is difficult in transmittance settings due to the poor contrast information obtained by non-optimized source-detector configurations. Therefore, image reconstructions from linear source locations suffer not only from fewer photon counts but also from limited sampling of the different organs to be imaged. Even though engineering optimization such as SVA analysis will suggest that linear sampling provide a better conditioned inverse problem, non-linear probe settings is of paramount importance to imaging organs such as the heart.

In conclusion, this research presents computational evidence that the combination of time-gated pMC and spatial a-priori significantly improve the accuracy of diffuse optical tomography in pre-clinical imaging where the diffusion theory fails [37]. The proposed method is easy to implement, flexible and can be easily extended to fluorescemce molecular imaging. Thanks to the concurent development of a multi-spectral time-gated imager [33], we hope to use this new reconstruction strategy to construct an in-vivo functional atlas of murine model. Such an atlas will be extremely valuable for the biomedical community.

Acknowledgements

The authors are grateful to M. Guven and V. Venugopal for insightful technical discussions. X Intes gratefully acknowledges the technical support of the Computational Center for Nanotechnology Innovations (CCNI) at RPI.

References and links

1. J. S. Lewis, S. Achilefu, J. R. Garbow, R. Laforest, and M. J. Welch, “Small animal imaging. current technology and perspectives for oncological imaging,” Eur. J. Cancer 38(16), 2173–2188 (2002). [CrossRef] [PubMed]

2. V. Koo, P. W. Hamilton, and K. Williamson, “Non-invasive in vivo imaging in small animal research,” Cell. Oncol. 28(4), 127–139 (2006). [PubMed]

3. A. H. Hielscher, “Optical tomographic imaging of small animals,” Curr. Opin. Biotechnol. 16(1), 79–88 (2005). [CrossRef] [PubMed]

4. F. Azar, and X. Intes, Translational Multimodality Optical Imaging, chap. Introduction to Clinical Optical Imaging, p. 1 (Artech House, Norwood, 2008).

5. S. Bjoern, S. V. Patwardhan, and J. P. Culver, “The influence of heterogeneous optical properties on fluorescence diffusion tomography of small animals,” in 2006 OSA/BOSD, AOIMP, TLA (2006).

6. L. Hervé, A. Koenig, A. Da Silva, M. Berger, J. Boutet, J. M. Dinten, P. Peltié, and P. Rizo, “Noncontact fluorescence diffuse optical tomography of heterogeneous media,” Appl. Opt. 46(22), 4896–4906 (2007). [CrossRef] [PubMed]

7. S. L. Jacques and B. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt. 13(4), 041302 (2008). [CrossRef] [PubMed]

8. X. Intes and B. Chance, “Non-PET functional imaging techniques: optical,” Radiol. Clin. North Am. 43(1), 221–234, xii (2005). [CrossRef] [PubMed]

9. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Ganudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). [CrossRef]

10. S. R. Arridge and W. R. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998). [CrossRef]

11. X. Intes and B. Chance, “Multi-frequency diffuse optical tomography,” J. Mod. Opt. 52(15), 2139–2159 (2005). [CrossRef]

12. U. Burcin, O. Birgul, R. Shafiiha, G. Gulsen, and O. Nalcioglu, “Diffuse optical tomographic reconstruction using multifrequency,” J. Biomed. Opt. 11, 054008 (2006). [CrossRef]

13. B. Chance, M. Cope, E. Gratton, N. Ramanujam, and B. Tromberg, “Phase measurement of light absorption and scatter in human tissue,” Rev. Sci. Instrum. 69(10), 3457–3481 (1998). [CrossRef]

14. A. Pineda, M. Schweiger, S. Arridge, and H. Barrett, “Information content of data types in time-domain optical tomography,” J. Opt. Soc. Am. A 23(12), 2989–2996 (2006). [CrossRef]

15. E. Alerstam, S. Andersson-Engels, and T. Svensson, “Improved accuracy in time-resolved diffuse reflectance spectroscopy,” Opt. Express 16(14), 10440–10447 (2008). [CrossRef] [PubMed]

16. X. Intes, S. Djeziri, Z. Ichalalene, N. Mincu, Y. Wang, P. St-Jean, F. Lesage, D. Hall, D. Boas, M. Polyzos, D. Fleiszer, and B. Mesurolle, “Time-domain optical mammography SoftScan: initial results,” Acad. Radiol. 12(8), 934–947 (2005). [CrossRef] [PubMed]

17. S. Lam, F. Lesage, and X. Intes, “Time Domain Fluorescent Diffuse Optical Tomography: analytical expressions,” Opt. Express 13(7), 2263–2275 (2005). [CrossRef] [PubMed]

18. S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008). [CrossRef] [PubMed]

19. K. M. Yoo, F. Liu, and R. R. Alfano, “When does the diffusion approximation fail to describe photon transport in random media?” Phys. Rev. Lett. 64(22), 2647–2650 (1990). [CrossRef] [PubMed]

20. J. C. Rasmussen, T. Pan, A. Joshi, T. Wareing, J. McGhee, and E. M. Sevick-Muraca, “Comparison of radiative transport, Monte Carlo, and Diffusion forward models for small animal optical tomography,” in 2007 IEEE ISBI, pp. 824–827 (2007).

21. Y. P. Kumar and R. M. Vasu, “Reconstruction of optical properties of low-scattering tissue using derivative estimated through perturbation Monte-Carlo method,” J. Biomed. Opt. 9(5), 1002–1012 (2004). [CrossRef] [PubMed]

22. J. Chen, V. Venugopal, and X. Intes, “Diffuse optical tomography with time-gated perturbation Monte Carlo method,” Proc. SPIE 7171, 717113–717119 (2009). [CrossRef]

23. Translational Multimodality Optical Imaging, F. Azar and X. Intes (eds.), Artech House, Norwood, 2008.

24. X. Intes, C. Maloux, M. Guven, B. Yazici, and B. Chance,“Diffuse optical tomography with physiological and spatial a priori constraints,” Phys. Med. Biol. 49(12N155CN), 163 (2004). [CrossRef]

25. S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo Model of Light Propagation in Tissue,” in SPIE Institute Series IS 5, (Society of Photo-Optical Instrumentation Engineers, Bellingham, Wash., 1989), pp. 102–111.

26. L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47(2), 131–146 (1995). [CrossRef]

27. D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [PubMed]

28. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001). [CrossRef]

29. A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. 44(11), 2082–2093 (2005). [CrossRef] [PubMed]

30. M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50(12), 2837–2858 (2005). [CrossRef] [PubMed]

31. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007). [CrossRef] [PubMed]

32. D. Stout, P. Chow, R. Silverman, R. M. Leahy, X. Lewis, S. Gambhir, and A. Chatziioannou, “Creating a whole body digital mouse atlas with PET, CT and cryosection images,” Mol. Imaging Biol. 4(4), S27 (2002).

33. V. Venugopal, J. Chen, and X. Intes, “Quantifying optical properties in small animals using MR-guided multi-spectral time-resolved imaging,” Proc. SPIE 7171, 717114–717118 (2009). [CrossRef]

34. W. F. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]

35. G. Alexandrakis, F. R. Rannou, and A. 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). [CrossRef] [PubMed]

36. A. E. Cerussi, A. J. Berger, F. Bevilacqua, N. Shah, D. Jakubowski, J. Butler, R. F. Holcombe, and B. J. Tromberg, “Sources of absorption and scattering contrast for near-infrared optical mammography,” Acad. Radiol. 8(3), 211–218 (2001). [CrossRef] [PubMed]

37. J. Chen, and X. Intes, “Time-resolved perturbation Monte Carlo for 3D optical imaging in small animals,” 34th Annual Northeast Bioengineering Conference, April 4th 2008.

38. E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. 13(041), 304 (2008). [CrossRef]

	Background	Heart	Stomach	Liver	Kidneys	Lungs
BV [μM]	0.06	0.15	0.1	0.20	0.056	0.1
StO₂ (%)	75	75	50	75	75	85
H₂O (%)	50	50	80	70	80	85

$μ_{a}$ (cm $^{- 1}$ )	Background	Heart	Stomach	Liver	Kidneys	Lungs
700nm	0.217	0.538	0.563	0.718	0.205	0.281
750nm	0.251	0.607	0.536	0.811	0.243	0.371
800nm	0.268	0.655	0.438	0.874	0.257	0.450
850nm	0.332	0.798	0.501	1.065	0.324	0.574
900nm	0.384	0.909	0.579	1.214	0.381	0.664
950nm	0.534	1.041	0.796	1.402	0.628	0.928

	Data	n-Lin.	Baye.	Backg.	Heart	Liver	Kidneys	Stomach	Lungs
BV	CW			17.2±35.7	−25.3±14.4	24.4±40.5	−32.0±20.6	21.2±21.5	3.1±27.5
	CW		×	12.4±3.4	−26.2±0.8	28.3±3.1	−36.0±1.6	20.1±3.3	−4.0±1.4
	TG			15.7±32.3	−20.5±15.2	22.7±38.2	−18.2±19.3	12.8±28.1	6.0±22.4
			×	11.9±5.9	−23.8±1.6	25.9±5.0	−18.6±2.8	14.6±4.7	3.9±3.0
		×	×	2.3±5.4	3.6±1.7	10.6±4.0	−4.2±3.0	3.0±4.7	13.5±3.7
StO₂	CW			3.6±22.4	−7.5±4.0	28.1±13.8	0.9±70.9	4.1±11.1	−8.4±6.0
	CW		×	−3.0±0.4	−8.8±0.1	14.8±0.7	−8.8±0.5	11.1±1.2	−5.7±0.2
	TG			2.9±16.5	−4.3±5.3	15.2±10.1	5.3±58.3	2.8±11.9	−10.5±5.0
			×	1.8±0.5	−4.2±0.3	9.80±1.0	−4.3±0.8	6.9±1.0	−7.6±0.3
		×	×	0.5±5.4	−1.1±1.7	10.7±4.0	−1.8±3.0	0.5±4.7	1.3±3.7
H₂O	CW			−48.8±16.7	−58.0±23.8	−64.2±18.2	−48.4±44.3	−50.1±21.4	−68.8±22.8
	CW		×	−7.0±4.8	−21.8±10.0	−32.8±7.3	−27.8±5.4	−19.1±12.1	−25.7±2.5
	TG			−49.3±14.3	−49.1±9.2	−60.8±10.2	−39.5±10.2	−44.1±16.0	−38.9±5.5
			×	−17.0±2.8	−22.4±3.2	−30.8±1.1	−35.0±8.2	−10.8±9.6	−27.6±0.3
		×	×	−20.9±3.6	−25.7±1.3	−33.8±6.6	−32.4±-5.0	−5.9±4.7	−13.7±6.7

	Background	Heart	Stomach	Liver	Kidneys	Lungs
BV [μM]	0.06	0.15	0.1	0.20	0.056	0.1
StO₂ (%)	75	75	50	75	75	85
H₂O (%)	50	50	80	70	80	85

$μ_{a}$ (cm $^{- 1}$ )	Background	Heart	Stomach	Liver	Kidneys	Lungs
700nm	0.217	0.538	0.563	0.718	0.205	0.281
750nm	0.251	0.607	0.536	0.811	0.243	0.371
800nm	0.268	0.655	0.438	0.874	0.257	0.450
850nm	0.332	0.798	0.501	1.065	0.324	0.574
900nm	0.384	0.909	0.579	1.214	0.381	0.664
950nm	0.534	1.041	0.796	1.402	0.628	0.928

Time-gated perturbation Monte Carlo for whole body functional imaging in small animals

Abstract

1. Introduction

2. Methodology

2.1 Time-resolved perturbation Monte Carlo

2.2 Functional imaging

2.3 Inverse problem

3. Simulation configuration

3.1 Geometric model and optode settings

3.1.1 Computational efficiency

3.1.2 Anatomically accurate simulations

3.2 Computational settings

3.3 Functional and optical parameters

4 Results

4.1 Computational feasibility

4.2 Time-gated reconstruction

5. Discussion

Acknowledgements

References and links

Cited By

Figures (4)

Tables (3)

Equations (13)

Optics Express