Home | About | Journals | Submit | Contact Us | Français |

**|**Biomed Opt Express**|**v.1(3); 2010 October 1**|**PMC3018054

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Equation of radiative transfer as a light propagation model
- 3. Two-step and multi-spectral PDE-constrained models
- 4. Results
- 5. Conclusions
- References and links

Authors

Related links

Biomed Opt Express. 2010 October 1; 1(3): 812–824.

Published online 2010 September 8. doi: 10.1364/BOE.1.000812

PMCID: PMC3018054

Hyun Keol Kim,^{1} Molly Flexman,^{1} Darrell J. Yamashiro,^{2} Jessica J. Kandel,^{3} and Andreas H. Hielscher^{1,}^{4,}^{5,}^{*}

Received 2010 June 14; Revised 2010 August 25; Accepted 2010 September 7.

Copyright ©2010 Optical Society
of America

This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported License, which permits download and redistribution, provided that the original work is properly cited. This license restricts the article from being modified or used commercially.

This article has been cited by other articles in PMC.

We introduce a transport-theory-based PDE-constrained multispectral model for direct imaging of the spatial distributions of chromophores concentrations in biological tissue. The method solves the forward problem (boundary radiance at each wavelength) and the inverse problem (spatial distribution of chromophores concentrations), in an *all-at-once* manner in the framework of a reduced Hessian sequential quadratic programming method. To illustrate the code’s performance, we present numerical and experimental studies involving tumor bearing mice. It is shown that the PDE-constrained multispectral method accelerates the reconstruction process by up to 15 times compared to unconstrained reconstruction algorithms and provides more accurate results as compared to the so-called two-step approach to multi-wavelength imaging.

Over the past decade near-infrared optical tomography (NIROT) has emerged as a viable imaging modality to monitor physiologically important parameters such as concentrations of oxyhemoglobin [HbO_{2}] and deoxyhemoglobin [Hb] in tissues [1–8]. This technique extracts tissue chromophore concentrations by exploiting the absorption data at multiple wavelengths, based on the fact that absorbers in tissue such as hemoglobin, water and lipid are linearly correlated with the tissue absorption. Initially researchers employed the so-called two-step method to obtain the concentrations of these absorbers and other chromophores. In this approach, the absorption coefficients are retrieved sequentially for each wavelength and then combined to obtain the concentrations through a direct matrix inversion or a least-squared method. However, it is well known [4,6,7] that this two-step approach sometimes produces unreasonable results (e.g. negative chromophore concentrations), and can show strong cross-talk between chromophores. To overcome this problem, researchers introduced a multispectral method in which chromophore concentrations are directly reconstructed by utilizing spectral constraints between chromophores into the reconstruction process. It has been shown that these direct methods provide improved accuracy and mitigate cross-talk. Li et al. [6] proposed a multispectral method based on the perturbation theory and demonstrated that it improved the concentration images of various chromophores. Corlu et al. [4,7] explored wavelength optimization for the continuous-wave multispectral optical tomography and validated their results through numerical simulations and with experimental data.

Both Li et al. and Corlu et al., however, based their work on the diffusion approximation (DA) to the more generally applicable equation of radiative transfer (ERT). It is well known that the diffusion approximation becomes less accurate when applied to imaging of small-tissue volumes and is further compromised if highly absorbing objects or fluid-filled regions, which contain, for example, cerebrospinal or synovial fluids, are considered [9]. The use of the full ERT alleviates these problems and provides accurate prediction of intensity measurements for imaging of small-tissue geometries. However, employing the ERT with the conventional two-step approach leads to prohibitively long computation times compared to the DA version, since the two-step model requires a large number of repeated forward simulations (based on ERT) that should be done for each wavelength during optimization. Thus it remains a challenging task to develop computationally efficient, accurate transport-theory-based multispectral chromophore reconstruction codes. This motivates our implementation of a fast ERT-based multispectral model for direct accurate imaging of chromophores concentrations in small-tissue volumes.

In this work we present a novel approach called a PDE-constrained multispectral imaging method that enables us to achieve two goals: first this approach accelerates the ERT-based image reconstruction process, and second it improves the chromophore reconstruction accuracy. To evaluate the new method, we provide numerical simulations and experimental results involving tumor bearing mice. Throughout these studies our new method is compared to the performance of the conventional unconstrained two-step method.

The frequency domain forward model for light propagation in tissue can be accurately modeled by the radiative transfer equation [10,11] as

$$[\Omega +{\mu}_{a}(\lambda )+{\mu}_{s}(\lambda )+i\frac{\omega}{c}]{u}_{\lambda}(\overrightarrow{r},\Omega ,\omega )=\frac{{\mu}_{s}^{\lambda}}{4\pi}{\displaystyle \underset{4\pi}{\int}{u}_{\lambda}(\overrightarrow{r},{\Omega}^{+},\omega )}\Phi ({\Omega}^{+},\Omega )d{\Omega}^{+}+{Q}_{\lambda}(\overrightarrow{r},\Omega ,\omega )$$

(1)

where *u _{λ}* is the spectral radiation intensity in unit [W/cm

We use an unstructured node-centered finite-volume approach [11] for the angular domain and a discrete ordinates formulation for the spatial domain, in order to obtain the discretized form of ** Au = b** where the sparse matrix

The tissue chromophore reconstruction is based on the linear correlation between the absorption coefficient and the chromophore concentration. The known mapping function of the *i*-th chromophore concentration C* _{i}* to the absorption coefficient μ

$${\mu}_{a}(\lambda )={\displaystyle \sum _{i=1}^{{N}_{c}}{\epsilon}_{i}(\lambda ){C}_{i}}$$

(2)

where *N _{c}* is the number of total chromophores that contribute to the absorption at wavelength λ. In the matrix form, Eq. (2) is given as [

Reconstruction of tissue chromophore concentrations can be performed in two different ways: a two-step method or a one-step multispectral method. Both methods will be described in some detail next.

The traditional two-step method deals with reconstructing the absorption data μ_{a}(λ) for each wavelength λ which corresponds to minimizing the objective function defined at wavelength λ given as

$$\varphi [{\mu}_{a}({\lambda}_{n}),{\mu}_{s}({\lambda}_{n})]=\frac{1}{2}{\displaystyle \sum _{d=1}^{ND}{\left|Q{u}_{d}^{\lambda}-{z}_{d}^{\lambda}\right|}^{2}}.$$

(3)

where *Q* is the measurement operator that projects the angular intensity distribution ${u}_{d}^{\lambda}$ onto the measurable quantity (i.e., radiation flux at the surface), ${z}_{d}^{\lambda}$ is the spectral measurement made at surfaces and *ND* is the total number of measurements. Then these images are converted into the chromophore concentration images. For example, consider four main absorbers (HbO_{2}, Hb, Water and Lipid) with four wavelengths measurements. In this case, we can decompose four absorption coefficient maps into into the chromophore images by using Eq. (2) as

$$\left[\begin{array}{l}{\text{HbO}}_{\text{2}}\\ \text{Hb}\\ \text{Water}\\ \text{Lipid}\end{array}\right]={\left[\begin{array}{cccc}{\epsilon}_{{\text{HbO}}_{\text{2}}}^{{\lambda}_{1}}& {\epsilon}_{\text{Hb}}^{{\lambda}_{1}}& {\epsilon}_{\text{Water}}^{{\lambda}_{1}}& {\epsilon}_{\text{Lipid}}^{{\lambda}_{1}}\\ {\epsilon}_{{\text{HbO}}_{\text{2}}}^{{\lambda}_{2}}& {\epsilon}_{\text{Hb}}^{{\lambda}_{2}}& {\epsilon}_{\text{Water}}^{{\lambda}_{2}}& {\epsilon}_{\text{Lipid}}^{{\lambda}_{2}}\\ {\epsilon}_{{\text{HbO}}_{\text{2}}}^{{\lambda}_{3}}& {\epsilon}_{\text{Hb}}^{{\lambda}_{3}}& {\epsilon}_{\text{Water}}^{{\lambda}_{3}}& {\epsilon}_{\text{Lipid}}^{{\lambda}_{3}}\\ {\epsilon}_{{\text{HbO}}_{\text{2}}}^{{\lambda}_{4}}& {\epsilon}_{\text{Hb}}^{{\lambda}_{4}}& {\epsilon}_{\text{Water}}^{{\lambda}_{4}}& {\epsilon}_{\text{Lipid}}^{{\lambda}_{4}}\end{array}\right]}^{-1}\left[\begin{array}{c}{\mu}_{a}({\lambda}_{1})\\ \begin{array}{l}{\mu}_{a}({\lambda}_{2})\\ {\mu}_{a}({\lambda}_{3})\end{array}\\ {\mu}_{a}({\lambda}_{4})\end{array}\right]$$

(4)

If more than four wavelengths are used to obtain four chromophores, i.e., *N _{c}* <

The one-step multispectral method retrieves the spatial distributions of chromphore concentrations directly by employing the multispectral measurements simultaneously into reconstruction. This can be achieved by minimizing the sum of the objective functions defined at each wavelength λ as:

$$f({C}_{i})={\displaystyle \sum _{n=1}^{{N}_{\lambda}}\varphi [{\mu}_{a}({\lambda}_{n}),{\mu}_{s}({\lambda}_{n})]}=\frac{1}{2}{\displaystyle \sum _{n=1}^{{N}_{\lambda}}{\displaystyle \sum _{d=1}^{ND}{\left|Q{u}_{d}^{{\lambda}_{n}}-{z}_{d}^{{\lambda}_{n}}\right|}^{2}}}$$

(5)

where *N _{λ}* is the number of measurement wavelengths and

We introduce here the PDE-constrained multispectral model. In order to understand the difference between the PDE-constrained multispectral and standard PDE-constrained models, we begin with a brief review of the standard PDE-constrained optimization.

Consider the standard optical tomographic problem defined at each wavelength as:

$$\text{min}f({\mu}_{\lambda};{u}_{\lambda})=\frac{1}{2}{\left|Q{u}_{\lambda}-{z}_{\lambda}^{obs}\right|}^{2}+\beta R\text{subjectto}c({\mu}_{\lambda};{u}_{\lambda})=A{u}_{\lambda}-b=0$$

(6)

where ** u** and

The most common approach for solving Eq. (6) is to treat the forward variable ** u** as a dependent variable of the inverse variable

$$\text{minimize}{f}^{}$$

(7)

which is often referred to as “unconstrained” because equality ** Au = b** no longer appears in Eq. (7), i.e.,

Another approach to solve Eq. (6), is to treat the forward variable ** u** and the inverse variable

$$L({\mu}_{\lambda},{u}_{\lambda};{\eta}_{\lambda})\frac{1}{2}{\left|Q{u}_{\lambda}-{z}_{\lambda}^{obs}\right|}^{2}+{\eta}_{\lambda}^{T}(A{u}_{\lambda}-b)$$

(8)

Here ** η** is called the vector of Lagrange multipliers. The simultaneous solutions of forward and inverse problems can be achieved at points satisfying the so-called first order Karush-Khun-Tucker (KKT) conditions [11,15] where the gradient of

However, the PDE-constrained approach of this kind [Eq. (8)] cannot be used directly with the one-step multispectral model. The reason is that the standard PDE-constrained optimization method [Eq. (8)] works only for a single wavelength-dependent constrained problem [Eq. (6)]. Thus it cannot treat all wavelengths data simultaneously within a PDE-constrained optimization framework. Therefore the standard PDE-constrained optimization approach has to be reformulated into a one-step multispectral model so that all multispectral data can be utilized simultaneously into reconstruction, which is the focus of this study.

In the following we provide a description of the formulation of PDE-constrained multispectral optimization. The basic idea behind this new method is to combine the multispectral model [Eq. (5)] and the PDE-constrained approach [Eq. (8)]. In other words, to be able to use data from all wavelengths we minimize the sum of the objective function defined at each wavelength λ [Eq. (5)] within the framework of PDE-constrained optimization [Eq. (8)]. This is achieved with the following formulations for PDE-constrained multispectral optimization:

$$\text{min}f(x;{u}_{\lambda})=\frac{1}{2}{\displaystyle \sum _{\lambda}{\left|Q{u}_{\lambda}-{z}^{obs}\right|}^{2}}+\beta R(x)\text{subjectto}c(x;{u}_{\lambda})=A{u}_{\lambda}-b=0$$

(9)

$$L(x,{u}_{\lambda};{\eta}_{\lambda})\frac{1}{2}{\displaystyle \sum _{\lambda}{\left|Q{u}_{\lambda}-{z}_{\lambda}^{obs}\right|}^{2}}+{\displaystyle \sum _{\lambda}{\eta}_{\lambda}^{T}(A{u}_{\lambda}-b)}$$

(10)

Here the inverse variable becomes a vector ** x** of all unknown chromophore concentrations, i.e.,

In the PDE-constrained multispectral model, the reduced Hessian SQP scheme solves the optimization problem [Eq. (9)] by minimizing the quadratic approximation of the Lagrangian function given by Eq. (10) subject to the linearization of the forward equation as follows:

$$\text{min}\Delta {x}^{kT}{g}_{r}^{k}+\frac{1}{2}\Delta {x}^{kT}{H}_{r}^{k}\Delta {x}^{k}\text{subjectto}{C}^{k}\Delta {p}^{k}+{(A{u}_{\lambda}-b)}^{k}=0$$

(11)

where Δ*p* = (Δ*x*, Δ*u _{λ}*)

$$\begin{array}{l}\Delta {x}^{k}=-{\left({H}_{r}^{k}\right)}^{-1}{g}_{r}^{k}\\ \Delta {u}_{\lambda}^{k}={A}^{-1}\left[{(A{u}_{\lambda})}^{k}-b-{\left(A{u}_{\lambda}\right)}_{x}^{kT}\Delta {x}^{k}\right]\end{array}$$

(12)

where the matrix-vector product ${({H}_{r}^{k})}^{-1}{g}_{r}^{k}$ is directly approximated by using the limited-memory BFGS updating formula [16]. Here the reduced gradient ${g}_{r}^{k}$ is given for the *j*-th chromophore concentration as:

$${g}_{r,j}^{k}={\displaystyle \sum _{\lambda}{\epsilon}_{j}^{\lambda}\left({u}_{\lambda}^{kT}{A}_{{x}_{j}}^{T}{\eta}_{\lambda}^{k}\right)+}{{x}_{j}}_{{f}^{k}}$$

(13)

The Lagrangian multiplier ${\eta}_{\lambda}^{k}$ at each wavelength λ can be obtained by solving the corresponding adjoint equation, given as:

$${A}^{kT}{\eta}_{\lambda}^{k}=-{Q}^{T}{(Q{u}_{\lambda}^{k}-{z}_{\lambda}^{obs})}^{*}$$

(14)

where the operator (.)* denotes the complex conjugate of the complex vector. The global convergence of the PDE-constrained multispectral algorithm is ensured by a line search on the following ${l}_{1}$ merit function that somehow balances the aims of reducing the objective function and satisfying the constraints, given as:

$${\mathrm{{\eta}_{k}}}_{(}$$

(15)

which is chosen here for its simplicity and low computational cost. The merit function [Eq. (15)] is successively monitored to ensure the global progress towards the optimal forward and inverse solutions, while a line search is performed to find a step length ${\alpha}^{k}$ that can provide a sufficient decrease in the merit function. As a consequence, the new iterate is given as:

$$\begin{array}{l}{x}^{k+1}={x}^{k}+{\alpha}^{k}\Delta x\\ {u}_{\lambda}^{k+1}={u}_{\lambda}^{k}+{\alpha}^{k}\Delta {u}_{\lambda}^{k}.\end{array}$$

(16)

The algorithm as described in Eqs. (11)–(16) is used for our numerical and experimental studies that will be presented next.

In this section, we apply the newly developed PDE-constrained multispectral method to directly recover the chromphore concentrations inside the medium. The performance of the algorithm is evaluated in terms of CPU times and accuracy. Without loss of generality, the scattering coefficient was assumed to be constant over the medium, given by μ_{s}^{′}(λ) = Aλ^{-b} (A = 45000, b = 1.3). In both the numerical and the experimental studies, we focus on the reconstruction of the oxy-hemoglobin (HbO_{2}) and deoxy-hemoglobin (Hb) concentrations since these are two major absorbers in tissue closely related to disease progression and physiological process. Water concentration was assumed constant over the medium.

First, we consider two-dimensional circular enclosure, as shown in Fig. 1
, filled with an inhomogeneous medium. Two objects with a diameter of 0.2cm are embedded in the background: one differs only in [HbO_{2}] and the one only in [Hb]. The corresponding chromophore concentrations are given in Table 1
for the background medium and the objects. The background medium has [HbO_{2}] of 40[μM] and [Hb] of 20[μM]. The object 1 has [HbO_{2}] of 80 while the object 2 has [Hb] of 40[μM]. These values are chosen here to mimic absorption coefficients of 0.1~0.2 (cm^{−1}), which correspond to the typical transport regime.

Circular enclosure (left) and absorption spectra of HbO_{2} (right): the test medium has two objects of perturbation: one in [HbO_{2}] and one in [Hb].

Eight sources are located around the enclosure surface at regular intervals. The sources are isotropic-emitting into the medium and intensity-modulated with a frequency of 600 MHz. The 64 detectors are equally spaced around the boundary surface. This yields a total of 512 source-detector pairs for reconstruction. For noise-free simulated data *z _{d}* in Eq. (5), we used a mesh of 11538 triangle elements with a S

Since we look at the reconstruction of [HbO_{2}] and [Hb], we only need two wavelengths to distinguish between these two chromophores. However, different sets of two wavelengths may affect the reconstruction accuracy. To illustrate the influence of the wavelengths set chosen, we consider two different sets of measurement wavelengths which are selected here from commonly available diode lasers: one set with (650nm, 830nm), the other set with (760nm, 830nm). Another reason for this set comes from the fact that, as shown in Fig. 1, the molar extinction coefficients of HbO_{2} and Hb are different from each other for chosen wavelengths, which is thus believed to be favorable to reconstructing HbO_{2} and Hb concentrations.

We use two metrics for the reconstruction quality: a correlation coefficient *ρ*_{c} and a deviations factor *ρ*_{d}, as defined in [17]. Accordingly, the larger correlation coefficient, close to 1, and smaller deviation factor, close to zero, represents the higher accuracy. Also the correlation coefficient may be used to indicate the degree of cross-talk that can be measured as the HbO_{2} (Hb) concentration divided by the Hb (HbO_{2}) concentration in Hb (HbO_{2}) position of the known object. Besides the effects of different wavelength sets, we also investigate the CPU times in both PDE-constrained and conventional multispectral methods, which is another important metric for evaluation of the performance of the proposed algorithm. All reconstructions are performed using two methods: the PDE-constrained one-step multispectral method and the unconstrained two-step method. Note that the unconstrained two-step method is based on the quasi-Newton (BFGS) updating scheme [16,17]. Hereafter this unconstrained two-step method will be denoted by simply the conventional method. The results are given in Table 2
and in Figs. 2
and and33
.

Convergence history of the PDE-constrained multispectral method with iteration: (a) forward error; (b) inverse error.

Images of reconstructed [HbO_{2}] (top row) and [Hb] (bottom row), obtained with the PDE-constrained multispectral and unconstrained two-step methods for two different wavelength sets.

As shown in Table 2, the PDE-constrained method leads to a significant saving in the CPU time in all cases. For the (650nm, 830nm) set, the PDE-constrained multispectral method reaches convergence in 14.5 min while the unconstrained two-step method takes about 168 min to converge at the same convergence criterion used in the PDE-constrained multispectral method. Therefore, the PDE-constrained multispectral method reduces the computation time by a factor of about 12. We observe a similar time saving in the other cases of a different wavelength set (see Table 2). The (760nm, 830nm) set data converges in 9.4 min using the PDE-constrained method and in 84 min using the unconstrained two-step method, respectively. Therefore, in this case, the PDE-constrained multispectral method accelerates the reconstruction process by a factor of about 10.

The main reason for this significant reduction in the CPU time can be explained by the fact that the PDE-constrained multispectral method does not require the exact solution of the forward problem at each iteration of optimization until it converges to the optimal solution, as mentioned earlier. Indeed we do not solve the forward problem for each wavelength, and instead we solve the linearized forward problem Eq. (12) which is much inexpensive to solve) for each wavelength. In other words, even a loose tolerance of 10^{−2} can be used to solve for the linearized forward problem. Since we use a GMRES iterative solver, the usage of a loose tolerance (of 10^{−2}) takes a much smaller number of iterations to a stopping criterion than the tight tolerance (of 10^{−10}) used in the unconstrained two-step method. As a consequence, this leads to a significant time savings through the overall optimization process.

Figure 2(a) and 2(b) show the convergence history of the PDE-constrained multispectral method with the number of iterations. As expected, it is observed that the PDE-constrained multispectral scheme solves the forward and inverse problems simultaneously through decreasing *all* forward [Fig. 2(a)] and inverse [Fig. 2(b)] errors *at once* in each of optimization iterations.

Figure 3 shows the reconstructed concentration images of two chromophores HbO_{2} and Hb obtained with two wavelength sets. Each row corresponds to the reconstructed images of HbO_{2} and Hb concentrations, respectively, and the columns display the reconstructed images for wavelength sets 1 and 2, for the PDE-constrained and unconstrained two-step methods, respectively. In terms of image quality, the two methods with sets 1 and 2 give different results. For set 1, it can be seen from Table 2 and Fig. 3 that the PDE-constrained method gives more accurate results than the conventional method. The correlation factor of the PDE-constrained method is 0.62 and 0.71 for [HbO_{2}] and [Hb], respectively. This is almost 20% and 40% better than when the two-step method is used (*ρ*_{c} = 0.54 and 0.56, respectively). For set 2, both methods give similar results for both [HbO_{2}] and [Hb]. When both wavelength sets are compared against the method, the two-step method works better with the 760-830 set rather than the 650-830 set. On the other hand, the PDE-constrained method outperforms the unconstrained two-step method for both wavelength sets. In particular, it gives best results for the 650-830 set.

Also, we found that the two-step unconstrained method produces larger deviation factors due to overestimation or underestimation in either HbO_{2} or Hb or both of the two (see Fig. 3). More evidences supporting this statement can be found from results obtained with the (760nm, 830nm) set. As compared to the results with the (650nm, 830nm) set, the two methods give lower accuracy both with respect to the correlation coefficient and the deviation factor as given in Table 2. Especially some cross-talk between two hemoglobin concentrations is observed in both of the two methods: the false Hb (HbO_{2}) perturbation is retrieved in the position of object 2 Eq. (1) where the medium has only inhomogeneity in HbO_{2} (Hb). In addition to cross-talk, the conventional method reveals larger artifacts and overestimation or underestimation in both [HbO_{2}] and [Hb]. This nature of the (760nm, 830nm) set may be explained by the concept of a condition number *κ*(*E*) of the molar extinction coefficient matrix *E* of HbO_{2} and Hb that represents a measure of how ‘well-posed’ a problem is: the larger *κ*(*E*) is, the more ‘ill-posed’ is the system. In other words, if the condition number is large, even a small error in the estimated absorption coefficient may cause a large error in the reconstructed hemoglobin concentration, and vice versa. In our cases, the condition numbers for the (650nm, 830nm) and (760nm, 830nm) sets are 1.76 and 3.17 respectively. Therefore it is evident that the (760nm, 830nm) set exhibits the nature of more ill-posed noise-sensitive problem, consequently making it more difficult to distinguish between the two unknowns, as compared to the (650nm, 830nm) set.

In addition to numerical studies, we have started to explore the code’s performance using experimental data obtained from tumor bearing mice. For this experiment 10^{6} cultured human Ewing sarcoma cells engineered to express luciferase (SK-NEP1-luc) were implanted intra-renally in NCR nude mice and allowed to grow until the tumor size reached about 1g as determined by weekly bioluminescence measurements. The experimental data was acquired with a continuous-wave digital optical tomography system that illuminates the target simultaneously with two wavelengths (λ = 760 nm and 830nm) modulated at 5 kHz and 7 kHz respectively. The system’s 16 sources and 32 detectors are configured in two rings, separated by 1.25cm, around a 3.175cm Delrin cylinder. Each ring contains 8 sources and 16 detectors arranged in a source-detector-detector-source configuration. The animal, anesthetized using isofluorane gas, was placed in the cylinder with the kidney and tumor located between the two rings of sources and detectors [see Fig. 4(a)
]. Then 1% Intralipid fluid was used to fill the remaining space in the cylinder to help reduce edge effects. The system uses silicon photodiodes to detect the transmitted and reflected photons passing through the cylinder. A digital-signal-processor (DSP) chip provides fast demodulation and control of the system giving an imaging frame rate of over 5 frames per second. More detailed information about the measurement system can be found elsewhere [18].

Optical imaging probe with mouse positioned in the cylinder (a) and signal changes observed over a 5 day period for the sum of all measurements (b).

To monitor the hemodynamic effects in the growing tumor, four sets of optical tomographic images (1000 frames of data - approximately three minutes of imaging) were obtained at days 0, 1, 3, and 5 after the tumor size was determined to have reached 1g. A 300-point subsection of the 1000 frames data with minimal motion was selected and averaged to remove any respiratory and other noise. After the data was acquired with the mouse in the cylinder a reference set of data was acquired for a homogeneous medium of 1% Intralipid.

Since the instrument used for this study does not provide absolute measurements due to the unknown calibration errors such as photon loss in fibers, we defined the following objective function [19]:

$$\varphi [{\text{HbO}}_{\text{2}}\text{,Hb}]=\frac{1}{2}{\displaystyle \sum _{\lambda}{\displaystyle \sum _{s}{\displaystyle \sum _{d}{\left|\frac{{M}_{\text{tar},s,d}^{\lambda}}{{M}_{\text{ref,}s,d}^{\lambda}}{P}_{\text{ref},s,d}^{\lambda}-{P}_{\text{tar},s,d}^{\lambda}\right|}^{2}}}}$$

(17)

where the subscripts *s* and *d* denote the numbers of sources and detectors. ${M}_{\text{tar},s,d}^{\lambda}$ and ${M}_{\text{ref,}s,d}^{\lambda}$ denote the spectral measurements at wavelength λ for the target medium of unknown optical properties and the reference medium of known optical properties, respectively. ${P}_{\text{ref},s,d}^{\lambda}$ and ${P}_{\text{ref},s,d}^{\lambda}$ are the corresponding forward predictions for the reference medium of known optical properties and the target medium of unknown optical properties, respectively. Here 1% Intralipid matching fluid is used as the reference medium.

Figure 4(b) shows the changes in measured signal intensity as it develops over the 5-day period when the experiments were performed. The value plotted as a function of time is the sum of measured optical signal for all source-detectors pairs, normalized by the sum of the day 0 signal. In general, we always saw a decline in signal intensity in almost all source detector pairs over time, which is caused by the continuous growth of the tumor, which is accompanied by a well-understood increased in vascularization (and hence blood volume) in the field of view [20,21].

Using data from all source-detector pairs we furthermore retrieved the three-dimensional spatial distributions of oxy- and deoxy-hemoglobin, (HbO_{2}) and (Hb) as well as their sum the total hemoglobin (THb). As with the numerical studies, we compare the performance of the two methods (PDE-constrained multispectral and conventional two-step methods) on the experimental data. We measure the CPU times of the two methods, and since the exact distributions of chromophores are not known for the mouse, we just provide a qualitative assessment of the results based on the expected biological change.

The reconstructed images of [THb], [HbO_{2}] and [Hb] are shown in Figs. 5
and and66
. Note that Figs. 5 and and66 show the results for the 2cmx3cmx3cm volume in the cylinder of height 8cm and diameter 3.175cm. We present here a drawing of the 3D volume contour of the reconstructed values above the given threshold, which can provide a better way for visual inspection of the time-trace results than the 2D cross-section maps. In terms of the CPU times, the PDE-constrained multispectral method took approximately 2 hrs to converge, while the conventional two-step method reached convergence in about 24 hours on a Dual Core Intel Xeon 3.33GHz processor: this constitutes acceleration factor of about 10. Furthermore, the two methods gave different trends of reconstructed chromophores over time. As can be seen in Fig. 5, the images generated by the PDE-constrained multispectral method clearly show that [HbO_{2}], [Hb], and [THb] increase in volume and value over time. This is expected as the tumor volume and vascular density increase over time as well [20,21]. This is also in agreement with the observed changes in optical signal shown in Fig. 4. On the other hand, this trend is not seen in the results obtained with the two-step method (Fig. 6). The [HbO_{2}] and [Hb] values decrease at 72 hrs and 5 days, respectively. We believe that this is an artifact introduced by the two-step method. Since the two-step method does not take any spectral constraints into account, it has been observed to produce less reliable results before [4,6,7] (e.g., artifacts, cross-talk and negative values).

Reconstructed distributions of [THb] (1st row), [HbO_{2}](2nd row) and [Hb] (3rd row) as observed over the 5 day period, using the proposed PDE-constrained multispectral method, with a wavelength set of 760nm and 830 nm. Shown is the 2cmx3cmx3cm volume that **...**

In this study, we present the first ERT-based PDE-constrained multispectral inverse method for direct imaging of chromophore concentrations in biological tissue. The proposed method is based on the reduced Hessian sequential quadratic programming method and solves the forward problem for each wavelength and the inverse problem for obtaining chromophore concentrations, in an *all-at-once* manner, which leads to a significant saving in the total image reconstruction time. We have evaluated the performance of the proposed algorithm through numerical studies and with experimental data obtained from tumor bearing mice.

We found that the PDE-constrained method not only reduces the image reconstruction time by a factor of up to 15 but also gives more accurate results as compared to the conventional two-step method. We also observed that the wavelengths pair (650 nm, 830 nm) produces more accurate results with less cross-talk than the wavelength pair (760 nm, 830 nm). The experimental results demonstrate that the PDE-constrained method gives reasonable results consistent with expected signal trends.

This work was supported in part by grants from the National Institutes of Health (National Cancer Institute grants 4R33CA118306 and U54CA126513-039001). Furthermore, Molly Flexman is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

1. L.C. Enfield, A.P. Gibson, J.C. Hebden, and M. Douek, “Optical tomography of breast cancer—monitoring response to primary medical therapy,” Targ. Oncol. (2009) at DOI 10.1007/s11523-009-0115-z.10.1007/s11523-009-0115-z [PubMed] [Cross Ref]

2. Boverman G., Fang Q., Carp S. A., Miller E. L., Brooks D. H., Selb J., Moore R. H., Kopans D. B., Boas D. A., “Spatio-temporal imaging of the hemoglobin in the compressed breast with diffuse optical tomography,” Phys. Med. Biol. 52(12), 3619–3641 (2007).10.1088/0031-9155/52/12/018 [PubMed] [Cross Ref]

3. J. Masciotti, F. Provenzano, J. Papa, A. Klose, J. Hur, X. Gu, D. Yamashiro, J. Kandel, and A. H. Hielscher, “Monitoring tumor growth and treatment in small animals with magnetic resonance optical tomographic imaging,” in Multimodal Biomedical Imaging; Fred S. Azar, Dimitris N. Metaxas, eds., SPIE International Symposium on Biomedical Optics, Proc. SPIE 6081, #608105 (2006).

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

5. Bluestone A. Y., Stewart M., Lasker J., Abdoulaev G. S., Hielscher A. H., Abdoulaev G. S., Hielscher A. H., “Three-dimensional optical tomographic brain imaging in small animals, part 1: hypercapnia,” J. Biomed. Opt. 9(5), 1046–1062 (2004).10.1117/1.1784471 [PubMed] [Cross Ref]

6. Li A., Zhang Q., Culver J. P., Miller E. L., Boas D. A., “Reconstructing chromosphere concentration images directly by continuous-wave diffuse optical tomography,” Opt. Lett. 29(3), 256–258 (2004).10.1364/OL.29.000256 [PubMed] [Cross Ref]

7. Corlu A., Durduran T., Choe R., Schweiger M., Hillman E. M. C., Arridge S. R., Yodh A. G., “Uniqueness and wavelength optimization in continuous-wave multispectral diffuse optical tomography,” Opt. Lett. 28(23), 2339–2341 (2003).10.1364/OL.28.002339 [PubMed] [Cross Ref]

8. Hielscher A. H., Bluestone A. Y., Abdoulaev G. S., Klose A. D., Lasker J., Stewart M., Netz U., Beuthan J., “Near-infrared diffuse optical tomography,” Dis. Markers 18(5-6), 313–337 (2002). [PubMed]

9. Hielscher A. H., Alcouffe R. E., Barbour R. L., “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43(5), 1285–1302 (1998).10.1088/0031-9155/43/5/017 [PubMed] [Cross Ref]

10. Ren K., Bal G., Hielscher A. H., “Frequency domain optical tomography based on the equation of radiative transfer,” SIAM J. Sci. Comput. 28(4), 1463–1489 (2006).10.1137/040619193 [Cross Ref]

11. Kim H. K., Hielscher A. H., “A PDE-constrained SQP algorithm for optical tomography based on the frequency-domain equation of radiative transfer,” Inverse Probl. 25(1), 015010 (2009).10.1088/0266-5611/25/1/015010 [Cross Ref]

12. M. Modest, Riative heat transfer, MacGraw-Hill Inc., New York, 1993.

13. Saad Y., Schultz M. H., “GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Comput. 7(3), 856–869 (1986).10.1137/0907058 [Cross Ref]

14. S. Prahl, “Optical properties spectra,” retrieved 16 March 2003, http://omlc.ogi.edu/spectra/index.html, 2001.

15. Abdoulaev G. S., Ren K., Hielscher A. H., “Optical tomography as a PDE-constrained optimization problem,” Inverse Probl. 21(5), 1507–1530 (2005).10.1088/0266-5611/21/5/002 [Cross Ref]

16. J. Nocedal and S. J. Wright, *Numerical optimization*, (Spinger-Verlag, New York, 1999).

17. Klose A. D., Hielscher A. H., “Quasi-newton methods in optical tomographic image reconstruction,” Inverse Probl. 19(2), 387 (2003).10.1088/0266-5611/19/2/309 [Cross Ref]

18. Lasker J. M., Masciotti J. M., Schoenecker M., Schmitz C. H., Hielscher A. H., “Digital-signal-processor-based dynamic imaging system for optical tomography,” Rev. Sci. Instrum. 78(8), 083706 (2007).10.1063/1.2769577 [PubMed] [Cross Ref]

19. Pei Y., Graber H. L., Barbour R. L., “Influence of systematic errors in reference states on image quality and on stability of derived information for DC optical imaging,” Appl. Opt. 40(31), 5755–5769 (2001).10.1364/AO.40.005755 [PubMed] [Cross Ref]

20. Yokoi A., McCrudden K. W., Huang J., Kim E. S., Soffer S. Z., Frischer J. S., Serur A., New T., Yuan J., Mansukhani M., O’toole K., Yamashiro D. J., Kandel J. J., “Human epidermal growth factor receptor signaling contributes to tumor growth via angiogenesis in her2/neu-expressing experimental Wilms’ tumor,” J. Pediatr. Surg. 38(11), 1569–1573 (2003).10.1016/S0022-3468(03)00562-1 [PubMed] [Cross Ref]

Articles from Biomedical Optics Express are provided here courtesy of **Optical Society of America**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |