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 (HbO2) 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.
4.1 Numerical study
First, we consider two-dimensional circular enclosure, as shown in
, filled with an inhomogeneous medium. Two objects with a diameter of 0.2cm are embedded in the background: one differs only in [HbO2] and the one only in [Hb]. The corresponding chromophore concentrations are given in
for the background medium and the objects. The background medium has [HbO2] of 40[μM] and [Hb] of 20[μM]. The object 1 has [HbO2] 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 HbO2 (right): the test medium has two objects of perturbation: one in [HbO2] and one in [Hb].
Chromophore concentrations and optical properties of the background medium and the objects
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 zd
in Eq. (5)
, we used a mesh of 11538 triangle elements with a S10
quadrature. The prediction Qud
in Eq. (5)
is obtained on a coarser mesh of 2854 triangular elements with a S6
quadrature. Measurements containing noise are simulated by adding an error term to zd
in the form
, where σ
is the standard deviation of measurement errors and
is the random variable with normal distribution. Here we chose a noise level of 15dB which represents the typical noise level encountered in optical tomography [17
Since we look at the reconstruction of [HbO2] 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 , the molar extinction coefficients of HbO2 and Hb are different from each other for chosen wavelengths, which is thus believed to be favorable to reconstructing HbO2 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 HbO2
(Hb) concentration divided by the Hb (HbO2
) concentration in Hb (HbO2
) 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
]. Hereafter this unconstrained two-step method will be denoted by simply the conventional method. The results are given in
Image quality of reconstructed [HbO2] and [Hb] obtained with the PDE-constrained multispectral and unconstrained two-step methods
Convergence history of the PDE-constrained multispectral method with iteration: (a) forward error; (b) inverse error.
Images of reconstructed [HbO2] (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 , 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 ). 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.
and 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  and inverse  errors at once in each of optimization iterations.
shows the reconstructed concentration images of two chromophores HbO2 and Hb obtained with two wavelength sets. Each row corresponds to the reconstructed images of HbO2 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 and 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 [HbO2] 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 [HbO2] 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 HbO2
or Hb or both of the two (see ). 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 . Especially some cross-talk between two hemoglobin concentrations is observed in both of the two methods: the false Hb (HbO2
) perturbation is retrieved in the position of object 2 Eq. (1)
where the medium has only inhomogeneity in HbO2
(Hb). In addition to cross-talk, the conventional method reveals larger artifacts and overestimation or underestimation in both [HbO2
] and [Hb]. This nature of the (760nm, 830nm) set may be explained by the concept of a condition number κ
) of the molar extinction coefficient matrix E
and Hb that represents a measure of how ‘well-posed’ a problem is: the larger κ
) 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.
4.2 Experimental study
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 106
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
]. 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
where the subscripts s
denote the numbers of sources and detectors.
denote the spectral measurements at wavelength λ for the target medium of unknown optical properties and the reference medium of known optical properties, respectively.
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.
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
Using data from all source-detector pairs we furthermore retrieved the three-dimensional spatial distributions of oxy- and deoxy-hemoglobin, (HbO2) 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], [HbO2
] and [Hb] are shown in
. Note that and 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 , the images generated by the PDE-constrained multispectral method clearly show that [HbO2
], [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
]. This is also in agreement with the observed changes in optical signal shown in . On the other hand, this trend is not seen in the results obtained with the two-step method (). The [HbO2
] 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
] (e.g., artifacts, cross-talk and negative values).
Fig. 5 Reconstructed distributions of [THb] (1st row), [HbO2](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 (more ...)
Fig. 6 Reconstructed distributions of [THb](1st row), [HbO2](2nd row) and [Hb](3rd row) over a 5-day period, using the conventional two-step method, with a wavelength set of 760nm and 830 nm. Shown is the 2cmx3cmx3cm volume that includes a growing tumor in vivo. (more ...)