Non-stationary reconstruction for dynamic fluorescence molecular tomography with extended kalman filter

Xin Liu; Hongkai Wang; Zhuangzhi Yan

doi:10.1364/BOE.7.004527

1. Introduction

Dynamic fluorescence molecular tomography (FMT) is important for bio-medical research since it allows non-invasive, tomographic imaging of the dynamic bio-distributions of fluorescent bio-markers within small animals in vivo [1–3]. With the development of optical contrast agents [4, 5], more applications can be expected in detecting drug efficacy, monitoring disease progression, and evaluating novel therapy, etc.

However, challenges remain in the reconstruction of dynamic FMT. The main reason is the current reconstruction methods [6–9] mainly focus on solving the stationary FMT tomographic problems. They are performed based on the (implicit) assumption that the fluorophore is stationary that it does not change during FMT imaging processes. Generally, the imaging time is needed to collect a complete set of measurements, which will be used as input data of the conventional reconstruction methods. However, in some applications, e.g., the drug delivery research, the fluorophore's concentration and location may change significantly when collecting these projection data. If the conventional reconstruction methods are directly applied to these time-varying fluorescence measurements, the reconstructed results may suffer from a high level of artifacts [10, 11]. In addition, based on these stationary methods, only one tomographic image can be obtained after scanning one circle projection data. Considering the rotation and the exposure time in in vivo experiments, the acquisition time is generally >1 min [3]. As a result, some fast changes of fluorophore in imaged object may not be detected during FMT imaging processes. To some extent, the above factors limit the widespread application of dynamic FMT in bio-medical research.

To address the problem, some theoretical and experimental solutions have been proposed [10–14]. Considering that during FMT imaging processes, the fluorophore may be assumed as a discrete time stochastic process, it is possible to solve the timing-varying fluorescence tomographic problems by the state-estimation approach [11–13]. As demonstrated by S. Prince and V. Kolehmainen et al. [11, 12], using the state-estimation method combined with extend kalman filter (EKF) technique [15], they successfully imaged the kinetic behaviors of fluorophore within imaged object. Despite of these successful applications, the current reconstruction method based on kalman filter was demonstrated by the fiber-based systems of slab geometry. Compared with the widely-used full-angle FMT imaging system, the fiber-based systems have some limitations in spatial sampling and experimental simplicity [16]. In addition, in the current EKF methods, the kalman gain was calculated based on the direct inverse matrix of the forward model of FMT. However, in FMT, the forward model matrix is highly ill-posed due to the high degree of light scattering in biological tissue. Together, the kalman gain is a core step in kalman filter processes. As a result, the imaging performance of FMT obtained by kalman filter may be affected if the kalman gain is calculated by the direct inverse method.

To overcome these limitations, in this paper, we apply kalman filter to time-varying fluorescence measurements which are acquired from a hybrid free-space, full-angle FMT/x-ray computed tomography (XCT) imaging system. Especially, to improve the imaging performance of FMT obtained by kalman filter, the generalized inverse of kalman gain is calculated by a second-order iterative method. The simulation, phantoms, and in vivo experiments are performed to evaluate the performance of the proposed method. Briefly, in the numerical simulation, two cases are demonstrated. Firstly, based on a three-dimensional (3-D) digital mouse model [17], we simulate the time-varying bio-distributions of indocyanine green (ICG) in the mouse's kidneys. In the second case, we simulate the movement of the fluorophore in the 3-D digital mouse during FMT imaging processes. For the phantom study, one tube filled with ICG is placed in one glass cylinder, which is used to imitate the timing-varying fluorescence imaging processes. In the in vivo experiment, utilizing the hybrid imaging system, we further image the bio-distributions of ICG in a nude mouse liver after the intravenous injection. The experimental results indicate that by using the proposed EKF-based second-order iterative (EKF-SOI) method, we can clearly resolve the dynamic bio-distributions of fluorophore within imaged object. In addition, the reconstruction time resolution can be greatly reduced (~2.5 sec/frame) in contrast to the conventional stationary reconstruction method, which makes it possible to detect the movable fluorophore during the imaging processes.

The outline of this paper is as follows. In section 2, the methods are detailed. In section 3, we describe the experimental materials. In section 4, the simulation, phantom, and in vivo experimental results are shown. Finally, we discuss the results and draw conclusions in section 5.

2. Methods

2.1. FMT imaging model

For the continuous wave FMT with point excitation source, the fluorescence measurement, $Φ$ , acquired from surface of the imaged object can be generated based on the first order Born approximation [18]

Φ (r_{d}, r_{s}) = Θ \int_{V} G_{}^{λ_{f l}} (r_{d}, r_{p}) ρ (r_{p}) G_{}^{λ_{e x c}} (r_{p}, r_{s}) d r_{p}

where

ρ (r_{p})

is the fluorescence yield to be reconstructed.

G_{}^{λ_{f l}} (r_{d}, r_{p})

and

G_{}^{λ_{e x c}} (r_{p}, r_{s})

is the Green's function describing the light propagating at the fluorescence (emission) and excitation wavelength, respectively.

Θ

is a calibration parameter which depends on generic system characteristics. With the known optical properties, these Green's functions can be obtained by solving the following diffusion equation using finite element method [9]

- \nabla \cdot [D (r) \nabla G (r)] + μ_{a} (r) G (r) = δ (r - r_{s})

where

μ_{a}

is the absorption coefficient,

D

is the diffusion coefficient, and

δ (r - r_{s})

is the point excitation source. By discretization, Eq. (1) can be rewritten in following form

Φ = ℘ (ρ)

where

℘ (\cdot)

is a forward model used to map the unknown fluorescence distribution to known measurements. In this paper, the diffusion equation is solved using COMSOL Multiphysics 3.3 (COMSOL, Inc., Burlington, MA, USA).

2.2. Reconstruction method based on EKF

For the conventional methods, in reconstruction processes, the fluorophore $ρ$ does not change during the acquisition of the measurement vector $Φ$ . However, in some applications, the fluorophore's concentration and location may change significantly during imaging processes, as mentioned above. If these conventional reconstruction stationary methods are used to solve the time-varying measurements, on one hand, the reconstructed results may suffer from a high level of artifacts [11, 12]. On the other hand, only one tomographic image can be obtained after scanning one circle projection data. As a result, it is impossible to image the movement of fluorophore during the imaging processes.

To address the problem, the state space approach is used in this paper. Briefly, in FMT imaging, the fluorophore $ρ$ can be assumed as a discrete time stochastic Markov process, as follows,

ρ_{t} = ℑ_{t - 1} ρ_{t - 1} + n_{t - 1}

where

ℑ_{t}

is a state transition matrix;

n_{t}

is a zero mean state noise process with a covariance matrix

Q_{t}

; t is a time steps (state index).

In addition, the fluorescence measurements can be derived from $ρ_{t}$ using the following forward model,

Φ_{t} = ℘_{t} (ρ_{t}) + m_{t}

where

Φ_{t}

is a measurement acquired at time t;

℘_{t} (\cdot)

is a forward model (observation matrix) corresponding to

ρ_{t}

;

m_{t}

is an observation noise with a covariance matrix

R_{t}

.

Equations (4) and (5) constitute the state-space representation of the time-varying FMT problem. In the linear case, the kalman filter, as a typical algorithm, can be used to solve the problem. However, in FMT, the forward model is generally a nonlinear function. For the nonlinear problem, the kalman filter is not applicable. To address this problem, the nonlinear Eq. (5) is firstly approximated to the following linear equation,

Φ_{t} \approx ℘_{t} (ρ_{t}^{*}) + J_{t} (ρ_{t}^{*}) (ρ_{t}^{} - ρ_{t}^{*}) + m_{t}

where

ρ_{t}^{*}

is a linearization point and

J = \partial ℘ (ρ) / \partial ρ

is a Jacobian matrix. After that, the kalman filter method is used. This method is called the extended kalman filter. The full equations for the (extend) kalman filter is given as follows,

State prediction : ρ_{t}^{(-)} = ℑ_{t - 1} ρ_{t - 1}

Variance prediction : C_{t}^{(-)} = ℑ_{t - 1} C_{t - 1} ℑ_{t - 1}^{T} + Q_{t - 1}

Kalman gain : G_{t} = C_{t}^{(-)} J_{t}^{T} {(J_{t} C_{t}^{(-)} J_{t}^{T} + R_{t})}^{- 1}

State update : ρ_{t} = ρ_{t}^{(-)} + G_{t} (Φ_{t} - ℘_{t} (ρ_{t}^{(-)}))

Variance update : C_{t} = (I - G_{t} J_{t}) C_{t}^{(-)}

where

G_{t}

is a kalman gain;

C_{t}

is a estimation error covariance matrix for the kalman filter state vector;

Q_{t}

is a covariance matrix associated with noise vector

n_{t}

;

R_{t}

is a covariance matrix of the observation errors

m_{t}

.

Based on Eqs. (7)-(11), we find that the fluorophore $ρ_{t}$ is calculated from current and past data. In order to incorporate the whole data in the sequence, the fixed-interval kalman smoother is used, which is achieved by performing a backward run with the following equation

ρ_{t - 1 | N} = ρ_{t - 1} + ℵ_{t - 1} (ρ_{t | N} - ρ_{t}^{(-)})

where

ρ_{t | N}

is the kalman smoother estimate and N is a total number of time points. The smoothing gain

ℵ_{t - 1}

is generated as follows,

ℵ_{t - 1} = C_{t - 1} G_{t - 1}^{T} C_{t}^{(-) - 1}

2.3. Second-order iteration for kalmain gain

Let $A = J C J_{}^{T} + R$ , if the generalized inverse of $A$ can be obtained, the kalmain gain can be calculated by Eq. (9). As mentioned above, in FMT, the forward model $J$ is highly ill-posed. As a result, the imaging performance of FMT obtained by EKF may be affected, if the kalman gain is generated based on the direct inverse matrix of $A$ . To more effectively calculate the generalized inverse of $A$ , in this paper, a second-order iterative method [19, 20] is used as follows,

S_{k + 1} = S_{k} (2 I - A S_{k})

where

S_{0}

is chosen as

a A^{T}

[19], with

a = 1 / λ_{\max}

.

λ_{\max}

is a maximum eigenvalue of

A \cdot A^{T}

, where T represents transposition operation. Referring to [19], the series will be convergent to the generalized inverse of

A

when

k \to \infty

. The detailed derivation can be found in [19].

3. Materials

3.1. Experimental setup

The numerical simulation, phantom, and in vivo experiments were performed based on a hybrid FMT/XCT imaging system [21]. As shown in Fig. 1, the imaged object was fixed on a rotating stage, which allowed rotation of the imaged object for full-angle projection acquisitions. Around the rotation stage, a full-angle, free-space FMT and micro-XCT system were constructed, respectively. The detailed information on the imaging system could be found in [21].

Fig. 1 The hybrid FMT/XCT imaging system.

Download Full Size | PDF

3.2. Numerical simulation

In the numerical simulation, the imaged object was generated from a digital mouse model [17], as shown in Fig. 2. The absorption ( $μ_{a}$ ) and scattering ( $μ_{s}^{'}$ ) coefficients from [22] were assigned to different organs of the mouse (liver, kidneys, stomach, spleen, and pancreas), respectively, to simulate photons propagation in heterogeneous tissues. The optical properties outside these organs were regarded as homogeneous.

Fig. 2 Schematic diagram of the simulation model. (a) The mouse 3-D geometry model used in simulation studies. The red points in (a) represent the excitation point source positions, located at the height of 0.9 cm and 1.7 cm, respectively. For each excitation location, the fluorescence is measured from the opposite side within $150^{o}$ field of view. (b) The coronal slice of the 3-D digital mouse at x = 0.35 cm. (c) The intensity curve of ICG in mouse kidneys, which is acquired from [24]. Inset in (c) shows that ICG intensities for 24 fluorescence projections corresponding to 24 angles from circle 1 (acquired during 0 min to 1 min). Different colors correspond to actual boundaries of different organs (red: livers; green: kidneys; stomach: cyan; pancreas: magenta; spleen: blue; bones: yellow; white: surface).

Download Full Size | PDF

For FMT imaging, the collimated light beam [see red points in Fig. 2(a)] was modeled as an isotropic point source, located at one mean free path of photon transportation beneath the surface, at the heights of 0.9 cm and 1.7 cm, respectively. In the imaging processes, $360^{o}$ full view imaging was performed with 24 projections in a $15^{o}$ step.

3.2.1. Case 1: simulated the metabolic processes of ICG in kidneys of the mouse

With the imaging model, we firstly simulated the time-varying metabolism processes of ICG in the mouse's kidneys. Briefly, according to Fig. 2(c), different ICG intensities were set to the kidneys at every time points during measurement data collection process. For example, when collecting the first circle of fluorescence projection data (acquired during 0-60s), 24 different ICG intensities were set to the kidneys at each measured time point ( $t_{1} - t_{24}$ ). After that, fluorescence projection images from different angles were generated from Eq. (1).

Here, ICG intensity curve was acquired from [24]. Based on our in vivo experiments [3, 14], the interval of adjacent projection was set to 2.5 sec.

3.2.2. Case 2: simulated the movement of fluorophore in the mouse

In the case 2, we simulated the movement of fluorophore in the mouse during the imaging time ( $t_{1} - t_{24}$ ) which is needed to collect a complete set of data (e.g., 24 projection images in the case) used as the input data of the conventional reconstruction methods. In detail, as shown in Fig. 3, in $t_{1} - t_{6}$ , the fluorophore (see the black cylinder in Fig. 3) was assumed to locate in the position 1. From $t_{7}$ to $t_{14}$ , the fluorophore moved to the position 2. Finally, the fluorophore changed to the position 3 between $t_{15}$ to $t_{24}$ .

Fig. 3 Simulating the movement of fluorophore in mouse during FMT imaging time ( $t_{1} - t_{24}$ ). In the imaging processes, we needed to collect a complete set of data (24 projection images), which will be used as the input data of the conventional reconstruction methods. Briefly, In $t_{1} - t_{6}$ , the fluorophore is assumed to locate the position 1. In $t_{7} - t_{14}$ , the fluorophore moves to the position 2. From $t_{15}$ to $t_{24}$ , the fluorophore changes to the position 3.

Download Full Size | PDF

Here, a cylinders (2 mm in diameter and 4 mm in height) were used as fluorescence targets. In the simulation processes, the optical properties and the imaging parameters are the same with those used in the case 1.

3.3. Phantom experiment

The phantom experiment was performed to further validate the performance of the proposed method. As shown in Fig. 4(a), the phantom was made of a glass cylinder (outer diameter of 3.0 cm) filled with a mixture of intralipid, water, and ink, with an absorption coefficient of $μ_{a} =0 .3 {cm}^{-1}$ and a reduced scattering coefficient of $μ_{s}^{'} = 10 .0 {cm}^{-1}$ . The optical parameters are similar to the bulk optical properties of mouse tissues. A transparent glass tube filled with 1 $μM$ ICG [outer diameter of 4.6 mm, see red point in Fig. 4(a)] were immerged in the cylinder phantom and used as the fluorescence target.

Fig. 4 Setup for phantom experiment. (a) The phantom used in the experiment, made of a glass cylinder (outer diameter of 3.0 cm) containing the matching fluid (intralipid, ink, and water), with $μ_{a} =0 .3 {cm}^{-1}$ and $μ_{s}^{'} = 10 .0 {cm}^{-1}$ . A glass tube (outer diameter of 4.6 mm and height of 7 mm) [see red regions in (a)] filled with 1 $μM$ ICG is immerged inside the phantom and used as fluorescence target. The back point in (a) describes the excitation source location. (b) ICG concentration curve which is used to imitate the time-varying measurements during FMT imaging processes. Here, the ICG concentration curve is the same as that used in Fig. 2(c).

Download Full Size | PDF

Measurement data were obtained from the hybrid FMT/micro-XCT imaging system. In FMT imaging process, 24 fluorescence images were collected at every $15^{o}$ . Subsequently, 72 white light images were acquired in $5^{o}$ steps to form a 3-D geometry of the imaged phantom. Finally, a steel anchor point, which could be imaged in both two imaging systems, was used to provide the relative height information for registration.

Referring to [23], in FMT imaging processes, the fluorescence measurement $Φ$ is approximately linearly with the fluorescence yield $ρ$ . So, for simplification, in the phantom experiment, during the acquisition of a circle of data, the time-varying measurements are generated based on the following equation

{\hat{Φ}}_{k} = Φ_{k} \times (1+ ξ_{k})

where

Φ_{k}

is the acquired fluorescent photon density of the

k t h

projection angle in the phantom experiment;

{\hat{Φ}}_{k}

is the synthetics fluorescent photon density of the

k t h

projection angle;

ξ_{k}

is the ICG concentration at the

k t h

time point (angle), which is acquired from the concentration curve in Fig. 4(b). After obtaining the full-angle synthetics measurements (24 angles in the case), these data are assembled to form the time-varying fluorescence measurement sequence, which are used to as the input data of the reconstruction.

3.4. In vivo experiment

In the in vivo imaging, an eight-week-old nude mouse, injected with 100 $μ L$ of 50 $μ g / m L$ ICG via tail vein, was suspended on the rotation stage and then continuously rotated 120 circles to monitor the metabolism processes of ICG in the liver. The acquisition time for collecting one circle of projection data was about 1 min. For each circle, 24 fluorescence images were acquired in $15^{o}$ steps. To correct for optical inhomogeneities, 24 excitation light images were subsequently acquired.

After optical imaging, a hepatobiliary contrast agent (Fenestra LC) was slowly intravenously injected to enhance the anatomical structure of livers in XCT imaging. 60 minutes after the injection of the contrast agent, the mouse was scanned using the micro-XCT to obtain the anatomical structures of organs. Finally, similar to the phantom experiments, a steel anchor point was used to provide the registration information. Throughout experiments, anesthesia was maintained with an isoflurane veterinary vaporizer.

4. Results

4.1. Simulation study

4.1.1. Reconstructed the metabolic processes of ICG through the kidneys of the mouse

Figure 5 shows the reconstructed tomographic images, obtained by the conventional stationary reconstruction method and the proposed reconstruction method based EKF, respectively. Here, the 24 projection data from the first circle were used as the input data of the reconstruction methods. When collecting one frame data, the fluorophore's concentrations vary significantly. Figures 5(a) and 5(d) show the reconstruction results obtained by the conventional stationary method. Figures 5(b) and 5(e) show the reconstruction results obtained by the EKF method with the kalman gain calculated by the second-order iterative method, termed as EKF-SOI. In contrast, Figs. 5(c) and 5(f) show the reconstruction results obtained by the EKF with the kalman gain solved by the direct inverse matrix method, termed as EKF-DIM. The reconstruction results [Figs. 5(b), 5(c), 5(e), and 5(f)] are imaged at 60s (corresponding to the 24th projection angle).

Fig. 5 Comparison of the reconstructed results obtained by the conventional stationary reconstruction method and the proposed reconstruction method based EKF. Here, the 24 projection data from the first circle were used as the input data of the reconstruction methods. When collecting one frame data, the fluorophore's concentrations vary significantly. (a) and (d) The reconstruction results obtained by the conventional stationary method. (b) and (e) The reconstruction result obtained by the proposed EKF-SOI method, which are imaged at 60s (corresponding to the 24th projection angle). (c) and (f) The reconstruction result obtained by the EKF-DIM method at 60s. Different colors correspond to actual boundaries of different organs (red: livers; green: kidneys; stomach: cyan; pancreas: magenta; spleen: blue; bones: yellow; white: surface).

Download Full Size | PDF

For reconstruction, the volume of interest was a $2 .0 cm \times 3. 0 cm \times 2 .6 cm$ 3-D region and sampled to $21 \times 31 \times 27$ voxels. Only 7,432 voxels inside the imaged object and 28,976 source-detector pairs were used in the reconstruction processes. For the stationary method, the reconstruction was performed by ART method. For the reconstruction based on EKF, the initialization values of the EKF were set as $ρ_{0} = 0$ , $C_{0} = 0.025 I$ , $Q_{0} = 0.01 I$ , and $R_{0} = 0.0001 I$ , respectively. $ℑ_{t}$ was chosen as identity matrix. For the EKF-SOI reconstruction, the number of the second-order iteration was fixed to 20. The above parameters were selected experimentally based on the results. For the linearization point $ρ_{t}^{*}$ , we used the current value of the kalman predictor estimate at each time point, i.e., $ρ_{t}^{*} = ρ_{t}^{(-)}$ .

Figure 6 shows the reconstructed FMT image sequence at all time points ( $t_{1} - t_{24}$ , with a time interval 2.5s) during the data acquisition time for collecting the first circle of projection data (0-60s). Furthermore, Fig. 6(b) shows the 3-D rendering of the reconstructed FMT image sequence. Here, the entire tomographic image sequence (24 frame) is obtained by the proposed EKF-SOI method.

Fig. 6 The reconstructed FMT images (total in 24 frames, with a time interval 2.5 s) based on one circle of projection data (24 projection data), which are imaged at different time points (2.5s, 5s, 7.5s, 30s, 45s, and 60s). (a) The reconstructed 2-D tomographic sequence at different time points. (b) The 3-D rendering of the reconstructed images. These tomographic images are obtained by the proposed EKF-SOI method. All images are displayed in the same range, respectively.

Download Full Size | PDF

The results indicate that it is difficult to resolve the fluorescence distributions in the kidneys based on the conventional stationary method [see Figs. 5(a) and 5(d)]. In contrast, when extend kalman filter is used in reconstruction processes, we can clearly resolve the ICG bio-distributions in the kidneys, as shown in Figs. 5(b), 5(c), 5(e), and 5(f). Furthermore, it can be seen that in this case, the reconstruction results from the proposed EKF-DIM method are similar with those obtained by the EKF-SOI method. On the other hand, we can also find that by using the conventional reconstruction method, only one tomographic image [see Figs. 5(a) and 5(d)] can be obtained based on the acquired 24 projection data from the first circle. For comparison, based on the same data set, 24 tomographic images ( $t_{1} - t_{24}$ ) corresponding to different angle can be obtained by the proposed EKF-SOI reconstruction method (see Fig. 6). That is, based on the EKF-SOI method, we can obtain the entire tomographic images at all time points ( $t_{1} - t_{24}$ , with a time interval 2.5s) during FMT imaging processes. In other words, by using the proposed EKF-SOI method, the reconstruction time resolution of FMT is reduced from 1 min (the stationary reconstruction method) to ~2.5 sec, which will be helpful for observing the fast biological activities.

4.1.2. Reconstructed the movement of fluorophore in the mouse

Figure 7 compares the tomographic images obtained by the stationary, EKF-SOI, and EKF-DIM methods, respectively. Especially, during the acquisition time (0-60s) that is needed to collect a complete set of measurements for traditional reconstruction method, the fluorophore moved three times within the mouse. The detailed information can be found in Section 3.2.2. Figure 7(a) shows the reconstruction result obtained by the stationary method. Figures 7(b)-7(e) and Figs. 7(f)-7(i) show the reconstruction result obtained by the EKF-SOI and EKF-DIM method, respectively. In the reconstruction processes, the initialization values of the EKF were set as $ρ_{0} = 0$ , $C_{0} = 0.025 I$ , $Q_{0} = 0.01 I$ , and $R_{0} = 0.0001 I$ and $ℑ_{t}$ was chosen as identity matrix. In the EKF-SOI reconstruction, the number of the second-order iteration was set as 20. These parameters were the same as those used in case 1 (Section 4.1.1).

Fig. 7 Comparison of the tomographic images obtained by the stationary, EKF-SOI, and EKF-DIM reconstruction method, respectively. Especially, during the acquisition time (0-60s) that is needed to collect a complete set of measurements for traditional reconstruction method, the fluorophore moved three times within the mouse. The detailed information can be found in Section 3.2.2. (a) The reconstruction result obtained by the stationary method. (b)-(e) The reconstruction result obtained by the proposed EKF-SOI method, which are imaged at 2.5s, 15s, 35s, and 60s, respectively. (f)-(i) The reconstruction result obtained by EKF-DIM method, imaged at the above time points. The white circles in (a) and black circles in (b)-(i) depict the actual position of fluorophore within the mouse.

Download Full Size | PDF

The results indicate that based on the acquired 24 projection data, only one tomographic image was generated by the stationary reconstruction method. As a result, it is impossible to resolve the different positions of fluorophore during imaging processes [see Fig. 7(a)]. In contrast, when the proposed EKF-SOI method is used in reconstruction processes, based on the same projection data, the entire 24 tomographic images at all time points (0-60 s, with a time interval 2.5 s) can be generated [see Figs. 7(b)-7(e)]. Based on these tomographic images at different time points ( $t_{1} - t_{24}$ ), we can clearly resolve the different position information of fluorophore during imaging processes [Position 1: Fig. 7(c); Position 2: Fig. 7(d); Position 3: Fig. 7(e)]. In addition, similar to the case 1, the results also indicate that there is no significant difference between the EKF-SOI and EKF-DIM results (see the 2nd and 3rd rows of Fig. 7).

Further, Fig. 8 shows the 3-D rendering of the reconstructed FMT image sequence, which is obtained by the proposed EKF-SOI method. Based on these reconstructed 3-D images ( $t = 2 .5s, 15s, 35s, and 60s$ ), we can resolve the moveable fluorophore. For example, Fig. 8(b) indicates that the fluorophore locates in the position 1. Figure 8(c) indicates that the fluorophore moves to the position 2. Similarly, Fig. 8(d) indicates that the fluorophore further moves to the position 3.

Fig. 8 The 3-D rendering of the reconstructed images at different time points (2.5s, 15s, 35s, and 60s). These reconstructed images are obtained by the proposed EKF-SOI reconstruction method based on one circle of projection data. All images are displayed in the same range.

Download Full Size | PDF

4.2. Phantom experiment

Figure 9 depicts the reconstruction images from the phantom experiment. Figure 9(a) shows the reconstruction results obtained by the stationary method, where one circle of projection data was used as the input data of the reconstruction. Figures 9(b) and 9(c) show the reconstruction results obtained by the EKF-SOI and EKF-DIM method, respectively, which were imaged at 17.5 s (corresponding to the 7th projection angle).

Fig. 9 The reconstruction results from the phantom obtained by the conventional stationary and EKF methods, where the time-varying measurements from 24 angles were used as the input data of the reconstruction. (a) The reconstruction results obtained by the conventional method based on circle of projection data (24 projection data). (b) The reconstruction result obtained by the proposed EKF-SOI method, which is imaged at imaged at 17.5s (corresponding to the 7th projection angle). (c) The reconstruction result obtained by the EKF-DIM method, imaged at 17.5s. The green curves depict the phantom boundary obtained by back-projecting the 72 white light images [25], which are used to validate the accuracy of registration.

Download Full Size | PDF

For the phantom reconstruction, the volume of interest was a $3 .0 cm \times 3. 0 cm \times 3. 3 cm$ 3-D region and sampled to $31 \times 31 \times 17$ voxels. 9,680 voxels inside the imaged object and 16,181 source-detector pairs were used for reconstruction. Similar to the numerical simulation, for the conventional stationary methods, the reconstruction was performed by ART method. In the EKF reconstruction, the parameters of the EKF were $ρ_{0} = 0$ , $C_{0} = 0.025 I$ , $Q_{0} = 0.01 I$ , and $R_{0} = 0.0001 I$ , respectively. For the EKF-SOI reconstruction, the second-order iteration was terminated after 20 iterations. Note that the time-varying measurements from 24 angles (one circle) were used as the input data of the reconstruction.

Figure 10 shows the reconstructed FMT image sequence through all time points ( $t_{1} - t_{24}$ , with a time interval 2.5 s) during FMT imaging process (0-60s). Furthermore, the second row of Fig. 10 depicts the time series of 3-D renderings of the reconstructed FMT image sequence. The entire tomographic image sequence (24 frame) is obtained by the proposed EKF-SOI reconstruction method.

Fig. 10 The reconstructed FMT image sequence (total in 24 frames) from the phantom experiments based on one circle of projection data (24 projection data). The first shows the reconstructed 2-D tomographic sequence at different time points (2.5s, 15s, 30s, 45s, and 60s). The second row shows the 3-D rendering of the reconstructed images at the above time points. All images are displayed in the same range, respectively. Note that by the proposed EKF-SOI method, we can obtain entire tomographic images through all time points during FMT imaging processes.

Download Full Size | PDF

The experimental results indicate that when using the conventional stationary method, the reconstructed image contain a high level of artifact [see Fig. 9(a)]. In contrast, when the proposed EKF-SOI reconstruction method is used, we can accurately resolve the ICG distribution in the phantom [see Fig. 9(b)]. In addition, it can also be found that compared with the EKF-SOI method, the reconstruction result from the EKF-DIM method contain more artifacts [see Fig. 9(c)]. On the other hand, similar to the numerical simulation, if using the conventional reconstruction method, only one tomographic image can be obtained based on the acquired 24 projection data. For comparison, based on the same data set, if using the EKF-SOI reconstruction method, we can acquire the entire 24 tomographic images corresponding to 1-24 angle during the imaging processes (see Fig. 10).

4.3. In vivo experiment

Figure 11 compares the reconstruction results from the in vivo experiment, obtained by the stationary, EKF-SOI, and EKF-DIM methods, respectively. Figure 11(a) shows the reconstruction image obtained by the stationary method, where the first circle of projection data was used as the input data of the reconstruction. Figure 11(b) shows the reconstruction results obtained by the proposed EKF-SOI method, which was imaged at 60s (corresponding to the 24th projection angle). For comparison, Fig. 11(c) shows the reconstruction results obtained by the EKF-DIM method at 60s.

Fig. 11 Comparisons of the reconstruction results in the in vivo experiment obtained by the stationary, EKF-SOI, and EKF-DIM methods, respectively. (a) The reconstruction results obtained by the conventional method based on circle of projection data (24 projection data). (b) and (c) The reconstruction results obtained by the EKF-SOI and EKF-DIM methods, which are imaged at imaged at 60 s (corresponding to the 24th projection angle).

Download Full Size | PDF

For the in vivo reconstruction, the volume of interest was a $2 .0 cm \times 2. 2 cm \times 3. 1 cm$ 3-D region and sampled to $21 \times 23 \times 32$ voxels. 5,132 voxels inside the imaged object and 5,117 source-detector pairs were used for reconstruction. In addition, to improve the reconstruction quality in the in vivo experiment, the reconstruction was performed by incorporating a heterogeneous forward model. Here, the forward model was constructed by assigning the midrange average optical properties to the heart, lungs, liver, and bones, which could provide an acceptable error of reconstructions in the mouse chest [26]. Similar to the simulation and phantom study, for the conventional stationary method, the reconstruction was performed by ART method. In the EKF reconstruction, we used $ρ_{0} = 0$ , $C_{0} = 0.025 I$ , $Q_{0} = 0.01 I$ , and $R_{0} = 0.0001 I$ , respectively. Furthermore, in the EKF-SOI reconstruction, the number of the second-order iteration was fixed to 5. These parameters were selected experimentally based on the results.

Furthermore, Fig. 12 depicts the reconstruction images that are imaged after injecting ICG. Figures 12(b)-12(d) show the reconstruction results obtained by the conventional methods, which were imaged at 0~1 min, 1~2 min, and 2~3 min after the injection of ICG, respectively. Figures 12(e)-12(i) show the reconstructed result obtained by the proposed EKF-SOI method based on the first circle projection data, which were imaged at 2.5s, 5.0s, 15s, 17.5s, and 30s after the injection of ICG. Similarly, Figs. 12(j)-12(m) show the EKF-SOI reconstruction results from the second circle, imaged at 62.5s, 65s, 92.5s, and 120s, respectively.

Fig. 12 The reconstruction results from the in vivo experiment. (a) The x-ray computed tomography image of liver region. (b)-(d) The FMT reconstruction results obtained by the conventional methods, which are imaged at 0~1 min, 1~2 min, and 2~3 min after the injection of ICG, respectively. (e)-(i) The FMT reconstructed results obtained by the proposed EKF-SOI method based on the first projection data, which are imaged at 2.5s, 5.0s, 15s, 17.5s, and 30s, respectively, after the injection of ICG. (j)-(m) The EKF-SOI reconstructed results from the second circle, which are imaged at 62.5s, 65s, 92.5s, and 120s, respectively.

Download Full Size | PDF

Finally, Fig. 13 depicts the time series of 3-D FMT renderings acquired after the injection of ICG, which are obtained by the proposed EKF-SOI method. The reconstruction images at different time points (2.5s, 5s, 45s, 60s, and 62.5s) are selected to show the metabolic processes of ICG in the mouse liver.

Fig. 13 The 3-D rendering of the reconstructed images at different time points (2.5s, 5s, 45s, 60s, and 62.5s) from the in vivo experiment. These reconstructed images are obtained by the proposed EKF-SOI method. All images are displayed in the same range.

Download Full Size | PDF

By comparing with the anatomical structure of the liver provided by the XCT image, the experimental results indicate that the proposed EKF-SOI method can provide a good reconstruction quality [see Fig. 11(b)]. However, there is significant difference between the reconstruction results obtained by the EKF-SOI and EKF-DIM methods [see Fig. 11(b) and 11(c)]. Here, the reconstruction performance of EKF-DIM is significantly affected. On the other hand, similarly to the simulation and phantom experiments, if using the EKF-SOI method, based on one circle of projection data, we can obtain the entire 24-frame tomographic images corresponding to different projection angles (time points) (see Fig. 12 and Fig. 13). In contrast, if using the stationary method, based on the same data (one circle of projection data), only one tomographic image can be obtained [see Fig. 11(a)].

5. Discussion and conclusion

Dynamic FMT imaging allows fast and non-invasively resolving the 3-D distribution of fluorescence bio-markers within small animal in vivo, which is helpful for bio-medical research. However, reconstruction of the dynamic FMT is a challenging problem. The main reason is that the current reconstruction methods are performed based on the implicit assumption that the fluorophore's distributions are stationary during FMT imaging processes. To overcome the limitation, the main aim of this paper is to solve the non-stationary (time-varying) fluorescence tomography problem, which is achieved by extended kalman filter combined with second-order iterative method.

It could be observed from the numerical simulation, phantom, and in vivo experimental results that when the EKF-SOI method was used in reconstruction processes, we could clearly resolve the dynamic distributions of ICG in the imaged object (Figs. 5 and 9) compared with the conventional stationary method. In addition, we could also find that by the proposed EKF-SOI method, we could obtain the entire tomographic images at all time points based on one circle of projection data (Figs. 6, 10, and 13). For comparison, when using the conventional reconstruction method, only one tomographic image could be obtained based on the same projection data. That was, by using the proposed method, the reconstruction time resolution of FMT was reduced from 1 min (the stationary reconstruction method) to ~2.5 sec. As a result, by the proposed EKF-SOI method, it was possible to resolve the movement of fluorescence bio-markers within imaged object (Figs. 7 and 8) during the course of collecting one circle of projection data.

Based on the above results, we believe that the proposed EKF-SOI reconstruction method provides an attractive technique for solving the dynamic FMT problem. Firstly, this technique allows to achieve the entire tomographic images at all time points during FMT imaging process, which is helpful for imaging fast biological activities. In addition, the reconstruction time resolution of FMT only depends on the rotation time of the rotation stage rather than the data acquisition time for scanning one circle projection data. As a result, the imaging time resolution can be further reduced when a fast rotation speed is used. Furthermore, this method is fairly flexible to extend to other optical imaging modalities, e.g., photoacoustic tomography (PAT) [27] and x-ray luminescence computed tomography (XLCT) [28].

It should be noted that the reconstructed images at each time point obtained by EKF-SOI are based on the measurements corresponding to one angle projection data. It seems not obtain the reconstruction images from all angles, which may lead to a compromised spatial resolution and highly ill-posed inverse problem due to the limited projections. In addition, the initialization value of the EKF parameters (e.g., the fluorophore $ρ$ , the estimation error covariance matrix $C_{t}$ , the noise covariance matrix $Q_{t}$ , and the covariance matrix of the observation errors $R_{t}$ ), are important, which may affect the performance of the proposed method. However, in this paper, the optimal EKF parameters are not determined. For example, in this work, the same value of $R_{0}$ is used for both the numerical simulation and real experiments. Considering that the measurement error of the real experiment is expected to be larger than the numerical simulation, a larger $R_{0}$ value might lead to better results for the EKF-SOI method in the real experiments. Moreover, the computation of the kalman smoother requires the all estimates from the forward run. Together, the forward matrix of FMT is generally large. As a result, the processing procedure of EKF-SOI is time-consuming, which is hard to meet the demands of dynamic imaging. Furthermore, to simplify the experiments, in the in vivo experiment, a thorough investigation of the reconstruction accuracy is not yet conducted. So, it is difficult to evaluate the performance of the proposed method in the in vivo experiment. To validate the reconstruction results from in vivo experiment, the multispectral epi-illumination cryoslicing imaging experiment may be performed [29]. These related work will be investigated in our future study.

In conclusion, we have applied the EKF-SOI method to FMT reconstruction to solve the time-varying FMT problem. The simulation, phantom, and in vivo experimental results indicate that based on the proposed method, we cannot only resolve the metabolic processes of fluorophore within the imaged object, but also detect the movement of fluorophore during the fast imaging processes, which will be helpful for resolving fast kinetic behaviors within small animals in vivo.

Funding

National Natural Science Foundation of China (NSFC) (81371604, 61571281, 6167124).

References and links

1. S. Patwardhan, S. Bloch, S. Achilefu, and J. Culver, “Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice,” Opt. Express 13(7), 2564–2577 (2005). [CrossRef] [PubMed]

2. X. Zhang, F. Liu, S. Zuo, J. Zhang, J. Bai, and J. Luo, “Fast reconstruction of fluorophore concentration variation based on the derivation of the diffusion equation,” J. Opt. Soc. Am. A 32(11), 1993–2001 (2015). [CrossRef] [PubMed]

3. X. Liu, X. Guo, F. Liu, Y. Zhang, H. Zhang, G. Hu, and J. Bai, “Imaging of indocyanine green perfusion inmouse liver with fuorescence diffuse optical tomography,” IEEE Trans. Biomed. Eng. 58(8), 2139–2143 (2011). [CrossRef]

4. X. Gao, Y. Cui, R. M. Levenson, L. W. K. Chung, and S. Nie, “In vivo cancer targeting and imaging with semiconductor quantum dots,” Nat. Biotechnol. 22(8), 969–976 (2004). [CrossRef] [PubMed]

5. M. C. Pierce, D. J. Javier, and R. Richards-Kortum, “Optical contrast agents and imaging systems for detection and diagnosis of cancer,” Int. J. Cancer 123(9), 1979–1990 (2008). [CrossRef] [PubMed]

6. K. O. Vasquez, C. Casavant, and J. D. Peterson, “Quantitative whole body biodistribution of fluorescent-labeled agents by non-invasive tomographic imaging,” PLoS One 6(6), e20594 (2011). [CrossRef] [PubMed]

7. A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging 24(10), 1377–1386 (2005). [CrossRef] [PubMed]

8. S. Patwardhan, S. Bloch, S. Achilefu, and J. Culver, “Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice,” Opt. Express 13(7), 2564–2577 (2005). [CrossRef] [PubMed]

9. A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004). [CrossRef] [PubMed]

10. M. J. Eppstein, D. E. Dougherty, T. L. Troy, and E. M. Sevick-Muraca, “Biomedical optical tomography using dynamic parameterization and Bayesian conditioning on photon migration measurements,” Appl. Opt. 38(10), 2138–2150 (1999). [CrossRef] [PubMed]

11. S. Prince, V. Kolehmainen, J. P. Kaipio, M. A. Franceschini, D. Boas, and S. R. Arridge, “Time-series estimation of biological factors in optical diffusion tomography,” Phys. Med. Biol. 48(11), 1491–1504 (2003). [CrossRef] [PubMed]

12. V. Kolehmainen, S. Prince, S. R. Arridge, and J. P. Kaipio, “State-estimation approach to the nonstationary optical tomography problem,” J. Opt. Soc. Am. A 20(5), 876–889 (2003). [CrossRef] [PubMed]

13. P. Hiltunen, S. Sarkka, I. Nissila, A. Lajunen, and J. Lampinen, “State space regularization in the nonstationary inverse problem for diffuse optical tomography,” Inverse Probl. 27(2), 025009 (2011). [CrossRef]

14. X. Liu, B. Zhang, J. Luo, and J. Bai, “4-D reconstruction for dynamic fluorescence diffuse optical tomography,” IEEE Trans. Med. Imaging 31(11), 2120–2132 (2012). [CrossRef] [PubMed]

15. M. Grewal and A. Andrews, Kalman Filtering: Theory and Practice Using MATLAB (New York: Wiley, 2001).

16. N. Deliolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360 degrees geometry projections,” Opt. Lett. 32(4), 382–384 (2007). [CrossRef] [PubMed]

17. 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]

18. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21(2), 158–160 (1996). [CrossRef] [PubMed]

19. X. Song, X. Xiong, and J. Bai, “A fast reconstruction algorithm for fluorescence optical diffusion tomography based on preiteration,” Int. J. Biomed. Imaging 2007, 23219 (2007). [CrossRef] [PubMed]

20. H. Wang, C. Wang, and W. Yin, “A pre-iteration method for the inverse problem in electrical impedance tomography,” IEEE Trans. Instrum. Meas. 53(4), 1093–1096 (2004). [CrossRef]

21. X. Guo, X. Liu, X. Wang, F. Tian, F. Liu, B. Zhang, G. Hu, and J. Bai, “A combined fluorescence and microcomputed tomography system for small animal imaging,” IEEE Trans. Biomed. Eng. 57(12), 2876–2883 (2010). [CrossRef] [PubMed]

22. 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]

23. X. Liu, F. Liu, and J. Bai, “A linear correction for principal component analysis of dynamic fluorescence diffuse optical tomography images,” IEEE Trans. Biomed. Eng. 58(6), 1602–1611 (2011). [CrossRef] [PubMed]

24. E. M. C. Hillman and A. Moore, “All-optical anatomical co-registration for molecular imaging of small animals using dynamic contrast,” Nat. Photonics 1(9), 526–530 (2007). [CrossRef] [PubMed]

25. D. Wang, X. Liu, Y. Chen, and J. Bai, “In-vivo fluorescence molecular tomography based on optimal small animal surface reconstruction,” Chin. Opt. Lett. 8(1), 82–85 (2010). [CrossRef]

26. D. Hyde, R. Schulz, D. Brooks, E. Miller, and V. Ntziachristos, “Performance dependence of hybrid x-ray computed tomography/fluorescence molecular tomography on the optical forward problem,” J. Opt. Soc. Am. A 26(4), 919–923 (2009). [CrossRef] [PubMed]

27. L. V. Wang and S. Hu, “Photoacoustic tomography: In vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef] [PubMed]

28. X. Liu, Q. Liao, and H. Wang, “In vivo x-ray luminescence tomographic imaging with single-view data,” Opt. Lett. 38(22), 4530–4533 (2013). [CrossRef] [PubMed]

29. A. Sarantopoulos, G. Themelis, and V. Ntziachristos, “Imaging the bio-distribution of fluorescent probes using multispectral epi-illumination cryoslicing imaging,” Mol. Imaging Biol. 13(5), 874–885 (2011). [CrossRef] [PubMed]

Non-stationary reconstruction for dynamic fluorescence molecular tomography with extended kalman filter

Abstract

1. Introduction

2. Methods

2.1. FMT imaging model

2.2. Reconstruction method based on EKF

2.3. Second-order iteration for kalmain gain

3. Materials

3.1. Experimental setup

3.2. Numerical simulation

3.2.1. Case 1: simulated the metabolic processes of ICG in kidneys of the mouse

3.2.2. Case 2: simulated the movement of fluorophore in the mouse

3.3. Phantom experiment

3.4. In vivo experiment

4. Results

4.1. Simulation study

4.1.1. Reconstructed the metabolic processes of ICG through the kidneys of the mouse

4.1.2. Reconstructed the movement of fluorophore in the mouse

4.2. Phantom experiment

4.3. In vivo experiment

5. Discussion and conclusion

Funding

References and links

Cited By

Figures (13)

Equations (15)

Biomedical Optics Express