Direct reconstruction in CT-analogous pharmacokinetic diffuse fluorescence tomography: two-dimensional simulative and experimental validations

Xin Wang; Yanqi Zhang; Limin Zhang; Jiao Li; Zhongxing Zhou; Huijuan Zhao; Feng Gao

doi:10.1117/1.JBO.21.4.046007

19 April 2016 Direct reconstruction in CT-analogous pharmacokinetic diffuse fluorescence tomography: two-dimensional simulative and experimental validations

Xin Wang, Yanqi Zhang, Limin Zhang, Jiao Li, Zhongxing Zhou, Huijuan Zhao, Feng Gao

Author Affiliations +

Journal of Biomedical Optics, Vol. 21, Issue 4, 046007 (April 2016). https://doi.org/10.1117/1.JBO.21.4.046007

Abstract

We present a generalized strategy for direct reconstruction in pharmacokinetic diffuse fluorescence tomography (DFT) with CT-analogous scanning mode, which can accomplish one-step reconstruction of the indocyanine-green pharmacokinetic-rate images within in vivo small animals by incorporating the compartmental kinetic model into an adaptive extended Kalman filtering scheme and using an instantaneous sampling dataset. This scheme, compared with the established indirect and direct methods, eliminates the interim error of the DFT inversion and relaxes the expensive requirement of the instrument for obtaining highly time-resolved date-sets of complete 360 deg projections. The scheme is validated by two-dimensional simulations for the two-compartment model and pilot phantom experiments for the one-compartment model, suggesting that the proposed method can estimate the compartmental concentrations and the pharmacokinetic-rates simultaneously with a fair quantitative and localization accuracy, and is well suitable for cost-effective and dense-sampling instrumentation based on the highly-sensitive photon counting technique.

1. Introduction

Though the utility of enhanced static imaging based on either targeting probes (cell-surface receptor, enzyme and reporter gene, and so on) or blood-pooling agents [indocyanine green (ICG)] is obvious with respect to cancer diagnosis and image-guided surgery, recent studies have observed significant intratumor spatial heterogeneity in expression of specific molecules and/or the metabolic process of nonspecific agents, as well as variations in expression patterns and agent permeability over time, suggesting that a regional biopsy at a fixed time point may not be representative of the complete disease.¹ In response, in vivo dynamic imaging approaches, as combined with advanced kinetic analyses, have been proposed to comprehensively examine the metabolic profiles of tumors, in addition to the cell (targeted protein) density or vascularization extent, for better identification and delineation of cancerous tissues from healthy ones. In particular, the heightened uptake of ICG, the only untargeted optical agent approved for human use, in many tumors compared with normal tissue, coined the “enhanced permeability and retention effect”, was originally observed back in the 1930s, and is now leveraged in many chemotherapies and disease characterizations.¹^–³ As an economical but effective way of overcoming the limitations of the conventional static diffuse fluorescence tomography (DFT), ICG-based pharmacokinetic DFT methods have been attempted by nature to provide enhanced contrast and specificity for tumor diagnosis, malignancy staging, and drug delivery assessment, and so on.⁴^–⁶

The current strategies for image reconstruction in the pharmacokinetic DFT can be categorized into indirect and direct methods: the former extracts images of the pharmacokinetic-rates by fitting the kinetic model to the time sequences of the fluorophore concentration first acquired with the conventional DFT,⁵^–⁹ while the latter directly reconstructs the pharmacokinetic images by incorporating a kinetic model to the inversion procedure.¹⁰^–¹⁴ For fully characterizing the ICG pharmacokinetic process, the indirect method normally requires the reconstruction of the interim ICG concentrations at not only a reasonably high temporal-resolution, but also over a large contrast variation. This presents a significant challenge to the state-of-the-art instrumentation in both the detecting sensitivity and data-acquisition mode. Although a charge-coupled-device-based system can achieve fast data-acquisition of one projection, sampling a set of the complete 360 deg projections at the required temporal-resolution (frame interval) and with enough signal-to-noise ratio (SNR) over the whole contrast range remains problematic and costly.⁶ Another adversity in the indirect method is the interim errors due to the severe ill-posedness of the DFT inversion, primarily shown as a saturating trend in the reconstruction quantitativeness with the contrast increasing.¹⁵^–¹⁹ The propagation of the errors into the analysis procedure will inevitably distort the time-course, and consequently deviate the estimation of the pharmacokinetic-rates.²⁰

Many approaches have been proposed to improve the performance of the indirect strategy: Liu et al.²¹ have proposed a four-dimensional (4-D) dynamic DFT scheme with data-sets of simultaneous $360 \deg - projections / frame$ and simplified the 4-D inversion with a Karhunen–Loève transformation. Zhang et al.²² reconstructed the fluorophore concentration variation instead of the absolute fluorophore concentration through a linear approximation with data-sets of instantaneous $360 \deg - projections / frame$ . Although the reconstruction could be achieved with an improved accuracy in the first method, the expensive requirement for fast data-acquisition remains. While for the second method, a reasonably small frame interval is requested to mitigate the deviation due to the temporally piecewise linear approximation, which also put forward a high demand for detection sensitivity which otherwise might lead to a bias in the estimation of subsequent kinetic parameters.

To overcome the inherent shortcomings of the indirect methods, the direct methods are currently under intensive investigation: A simple way is to incorporate a biexponential kinetic model into the DFT inversion to express a kinetics-to-measurement map, and directly extract the pharmacokinetics-related decaying constants from the data-sets of instantaneous $360 \deg - projections / frame$ .¹⁰^–¹³ Since the accurate pharmacokinetic fundamental is lacking, the biexponential model is empirical and thus, has no direct relationship with the tissue physiology. In addition, the reconstruction is nonreal-time and computationally intensive due to simultaneous exploitation of the total data-sets. A more complex method is based on the two-compartment model for the fluorophore kinetics to construct the kinetics-to-measurement map, and successively recovers the pharmacokinetic-rates in the conventional extended Kalman filtering (EKF) update from the frames of simultaneous complete 360 deg-projections.¹⁴ Although the use of the two-compartment model may lead to physiologically-consistent reconstructions, a complete set of the simultaneous 360 deg projections is required at each updating stage of the conventional-EKF, implying a high request for the temporal-resolution of instrumentation, and in addition, the performance of the conventional-EKF adversely hinges on the inaccessible a priori initial state of the model, as described in our previous work.²⁰

In this paper, we propose a direct strategy for the pharmacokinetic DFT by using the adaptive-EKF scheme for the generalized multicompartment model with the instantaneous sampling within the generalized framework of a discrete-detector-based CT-scanning mode. This strategy further relaxes the restrictions on acquiring a simultaneously sampled projection and simplifies the system configuration. The proposed method is first validated by simulated data for the two-compartment model with different pharmacokinetic contrasts in a two-dimensional (2-D) circular geometry, and then by pilot phantom experiments for the one-compartment model using a specifically-developed four-channel photomultiplier tube (PMT)-based photon-counting DFT system with CT-scanning mode. Both the simulative and experimental results suggest that the proposed method has a high feasibility to achieve simultaneous reconstructions of the compartmental concentrations and the pharmacokinetic-rates with reasonable accuracy.

2. Methods

Figure 1 shows the generalized schematic of a CT-analogous pharmacokinetic DFT, where the $D$ detecting positions $r_{d}^{(s)}$ ( $d = 1,2, \dots, D$ ) are evenly distributed on the surface of the phantom from 101.25 deg to 258.75 deg, opposite to an incidence (0 deg) at the $S$ source positions $r_{s}$ ( $s = 1,2, \dots, S$ ), at a given imaging plane, and $N_{d}$ discrete detectors are employed for the projection acquisition in a multiplexing way at a sampling period of $Δ T$ .²³ With such a configuration, $L (= D / N_{d})$ simultaneous samplings are performed to collect one complete projection, i.e., the data at all the detecting positions with regard to the illumination of one source, with the assumption that $N_{d}$ is smaller than or equal to $D$ . Normally, $P$ sets of the complete 360 deg projections are required to sufficiently cover a kinetic range of $P \times S \times L \times Δ T$ . The measurements are used in the normalized Born formulation of the excitation flux, $I_{x} (r_{d}, r_{s}, t)$ , and emission flux, $I_{m} (r_{d}, r_{s}, t)$ , as follows:²⁴

Eq. (1)

I_{n b} (p, s, l) = [\begin{array}{l} \frac{I_{m} {r_{(l - 1) N_{d} + 1}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}}{I_{x} {r_{(l - 1) N_{d} + 1}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}} \\ \frac{I_{m} {r_{(l - 1) N_{d} + 2}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}}{I_{x} {r_{(l - 1) N_{d} + 2}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}} \\ ⋮ \\ \frac{I_{m} {r_{l N_{d}}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}}{I_{x} {r_{l N_{d}}^{(s)}, r_{s}, [(p - 1) S L + (s - 1) L + l] Δ T}} \end{array}],

where

p = 1,2, \dots, P

,

s = 1,2, \dots, S

, and

l = 1,2, \dots L

denotes the number of the complete 360 deg data-sets, source positions and samplings, respectively.

Fig. 1

CT-analog pharmacokinetic DFT: (a) sketch of phantom and the configuration of the source and detectors and (b) measurement process of CT-analog pharmacokinetic DFT.

2.1.

Kinetics-to-Measurement Mapping

A crucial step toward the direct extraction of the pharmacokinetic-rates is to explicitly associate the time-dependent samplings with the kinetic model. In this work, a more generalized multicompartment model is methodologically adopted that describes the mutual permeability of administrated ICG agent among the conceptually divided $N_{c}$ compartments as well as its drainage out of tissue. The generalized model describes an increasingly function-detailed pharmacokinetic process: e.g., in the four-compartment model, the tissue is assumed to be composed of capillary region, interstitial fluid region, parenchymal cell region, and intracellular binding site, whereas, in the widely used two-compartment model, the interstitial fluid region, parenchymal cell region, and intracellular binding site are combined into a single compartment referred to as extra-cellular extra-vascular space (EES), and in the simplest one-compartment model, all the four compartments are merged and only the ICG drainage out of tissue is considered.⁷

For the generalized multicompartment model, a discrete nonlinear state-space (kinetic) model is expressed as follows, with regard to both the compartmental concentrations and the kinetics-relevant parameters:²⁵^–²⁸

Eq. (2)

[\begin{array}{l} C (r, k + 1) \\ θ (r, k + 1) \end{array}] = [\begin{array}{l} K [θ (r, k)] C (r, k) \\ θ (r, k) \end{array}] + [\begin{array}{l} ω (r, k) \\ ς (r, k) \end{array}],

where

C (r, k) = {[C_{1} (r, k), C_{2} (r, k), \dots, C_{N_{c}} (r, k)]}^{T}

with

C_{i} (r, k)

being the ICG tissue concentration in the

i

’th compartment (

i = 1,2, \dots, N_{c}

) at time

k Δ T

(

k = 1,2, \dots, K

);

K [θ (r, k)]

is the system matrix that depends on the pharmacokinetic-rates in an exponential fashion

Eq. (3)

K (θ (r, k)) = [\begin{matrix} τ_{11} (r) & τ_{12} (r) & \dots & τ_{1 N_{c}} (r) \\ τ_{21} (r) & τ_{22} (r) & \dots & τ_{2 N_{c}} (r) \\ ⋮ & ⋮ \\ τ_{N_{c} 1} (r) & τ_{N_{c} 2} (r) & \dots & τ_{N_{c} N_{c}} (r) \end{matrix}] = \exp ([\begin{matrix} - \sum_{i \neq 1} K_{1 i} - K_{1} & K_{21} & \dots & K_{N_{c} 1} \\ K_{12} & - \sum_{i \neq 2} K_{2 i} - K_{2} & \dots & K_{N_{c} 2} \\ ⋮ & ⋮ \\ K_{1 N_{c}} & K_{2 N_{c}} & \dots & - \sum_{i \neq N_{c}} K_{N_{c} i} - K_{N_{c}} \end{matrix}] Δ T);

ω (r, k)

and

ς (r, k)

are independent zero-mean Gaussian white noise processes;

θ (r, k) = {[τ_{11} (r, k), τ_{12} (r, k), \dots, τ_{N_{c} N_{c}} (r, k)]}^{T}

denotes a set of interim kinetic parameters that directly relates to the pharmacokinetic-rates:

K_{i j}

characterizes the permeability of ICG from the

i

’th compartment to the

j

’th compartment and

K_{i}

describes the ICG drainage out of the tissue.

To obtain the kinetics-to-measurement mapping, we first define a state vector, $X = {[C^{T} (r_{1}), θ^{T} (r_{1}), \dots, C^{T} (r_{N}), θ^{T} (r_{N})]}^{T}$ , to merge the compartmental concentrations and the interim kinetic parameters at all the nodes $r_{n}$ ( $n = 1,2, \dots, N$ ), and then combine the spatially-resolved multicompartment model with the DFT imaging equation that maps the sequences of the concentration image to the boundary samplings. This finally results in a new kinetic model regarding $X$ ¹⁴^,²⁰^,²⁵^–²⁸

Eq. (4)

{\begin{cases} X (k + 1) = A (X (k)) X (k) + δ (k) \\ I_{n b} (k) = Λ (k) X (k) + γ (k) \end{cases},

where

A = diag [K (r_{1}), I, K (r_{2}), I, \dots, K (r_{N}), I]

, with

I

being an identity matrix;

Λ (k) = ϵ η \ln (10) W (k) E

is referred to as the kinetics-to-measurement mapping matrix, with

η (= 0.016)

and

ϵ (= 0.013 {mm}^{- 1} μ M^{- 1})

being the quantum efficiency and extinction coefficients of ICG,⁷^,¹⁴^,²⁰ respectively,

E

is the

N \times (N_{c} + N_{c}^{2}) N

matrix transforming

X

to the total concentration image

C_{T} (= C_{1} + C_{2} + \dots + C_{N_{c}})

with

C_{i} = {[C_{i} (r_{1}), \dots, C_{i} (r_{N})]}^{T}

, and

W (k)

the

N_{d} \times N

-dimensional weighting matrix of the DFT imaging equation;¹⁹

δ (k)

and

γ (k)

are again approximated as independent zero-mean Gaussian white noise processes with the covariance matrices

Q

and

R

, respectively.

2.2.

Adaptive-Extended Kalman Filtering

The adaptive-EKF has been effectively used in our previous work for the parameter estimation of the two-compartment based nonlinear kinetic model using simultaneous sampling mode of the complete 360 deg projections,²⁰ and its generalization to the multicompartment kinetic model in Eq. (2) using the new instantaneous sampling mode is direct, as described in Eq. (4). Table 1 summarizes the general framework of an adaptive-EKF for the multicompartment model, which essentially describes a recursive procedure of “prediction-gain-update”:¹⁴^,²⁰ First, one-step ahead predictions of the state vector and the error-covariance matrix, $\hat{X} (k | k - 1)$ and $P (k | k - 1)$ , are made according to their previous step $\hat{X} (k - 1)$ and $P (k - 1)$ . Then the Kalman gain and the innovation sequence at the current step, $G (k)$ and $d (k)$ , are determined, respectively, based on the two predictions, $\hat{X} (k | k - 1)$ and $P (k | k - 1)$ , as well as the currently calculated the kinetics-to-measurement mapping matrix $Λ (k)$ . Finally, the state vector $\hat{X} (k)$ , the error-covariance matrix $P (k)$ and the Jacobian matrix $J (k)$ are updated for the current step, respectively, through the Kalman gain $G (k)$ as well as a temporally varying forgetting-factor $λ (k)$ that compensates the inaccuracy of the initial states and emphasizes the effect of the current data by averaging the previous innovations inside a moving window.²⁰^,²⁹^,³⁰

Table 1

Adaptive-EKF algorithm.

①Prediction:

{\begin{cases} \hat{X} (k | k - 1) = A [\hat{X} (k - 1)] \hat{X} (k - 1) \\ P (k | k - 1) = λ (k) J (k - 1) P (k - 1) J^{T} (k - 1) + Q \end{cases}

②Kalman gain:

{\begin{cases} G (k) = P (k | k - 1) Λ^{T} (k) {[Λ (k) P (k | k - 1) Λ^{T} (k) + R]}^{- 1} \\ d (k) = I_{n b} (k) - Λ (k) \hat{X} (k | k - 1) \end{cases}

③Update:

{\begin{cases} \hat{X} (k) = \hat{X} (k | k - 1) + d (k) G (k) \\ P (k) = [I - G (k) Λ (k)] P (k | k - 1) \\ J (k) = \frac{\partial [A (X (k)) X (k)]}{\partial X (k)} \end{cases}

3. Validations

To evaluate the performance of the proposed scheme, we conducted both the simulations and phantom experiments, with the former validating the theoretical feasibility while the latter demonstrated the technological effectiveness, of the method.

3.1.

Simulative Investigations

In the simulations, the two-compartment model ( $N_{c} = 2$ ) is used that has been widely-applied in the previous studies.⁷^,¹⁴^,²⁰ Tissues in this model are assumed to be composed of plasma and the EES, and the system matrix in Eq. (3) thus reduces to

Eq. (5)

K (θ (r)) = [\begin{matrix} τ_{11} (r) & τ_{12} (r) \\ τ_{21} (r) & τ_{22} (r) \end{matrix}] = \exp ([\begin{matrix} - K_{p e} (r) - K_{p} & K_{e p} (r) \\ K_{p e} (r) & - K_{e p} (r) \end{matrix}] Δ T),

where

K_{p e}

,

K_{e p}

, and

K_{p}

are the pharmacokinetic-rates describing the ICG leakage into the EES, the ICG drainage out of the EES, and the ICG drainage out of tissue, respectively. Correspondingly, we have

C (r, k) = {[C_{p} (r, k), C_{e} (r, k)]}^{T}

and

θ (r) = {[τ_{11} (r), τ_{12} (r), τ_{21} (r), τ_{22} (r)]}^{T}

.²⁰

C_{p} (r, k)

and

C_{e} (r, k)

denote the ICG concentration in the plasma and the EES, respectively. Note that the heterogeneity of

K_{p e}

and

K_{e p}

is the same as the physical tumor location and shape, while

K_{p}

is reasonably assumed to be constant throughout the domain.¹⁴^,²⁰

A circular tissue-emulating geometry with a 30-mm diameter is used with a circular tumor-embedded region of 6-mm diameter, centered at ( $x = 5 mm$ , $y = 0 mm$ ), as shown in Fig. 2. The absorption and reduced scattering coefficients are set to $μ_{a} = 0.035 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , respectively, for both the excitation and emission wavelengths, i.e., 780 and 830 nm, respectively. The finite-element method (FEM) based diffusion approximation (DA) solution in 2-D is used for both the forward and inverse calculations, with a FEM mesh of 721 nodes and 1350 triangular elements.¹⁹

Fig. 2

Geometric sketch of the phantom used in the simulative validations.

For the data generation, the 720-s time-course of the ICG compartmental concentrations $C$ is first calculated according to the two-compartment model at a sampling period of $Δ T = 2.5 s$ and corrupted by the Gaussian white noise with a SNR of $α_{ω} = 40 dB$ , i.e., $ω (r) = K [θ (r)] C (r) 10^{- α_{ω} / 20} R_{N}$ ( $R_{N}$ is the Gaussian random number) in Eq. (3), based on the pharmacokinetic parameters in Table 2. The ICG is assumed all in the plasma with a concentration of $1.0 μ M$ at the starting time of $t = 0 s$ , i.e., $C (r_{n}, 0) = {[0,1]}^{T}$ . Then the outward flux on the boundary are simulated using the FEM-DA-based forward calculation for a reference illumination-detection configuration of $S = 16$ , $D = 8$ , and $N_{d} = 4$ , which is in agreement with the following experimental setup, and added a Gaussian noise with the SNRs of $α_{x} = 50 dB$ and $α_{m} = 40 dB$ for the excitation and emission, respectively.¹⁵^–¹⁹

Table 2

Pharmacokinetic-rates with three contrasts for simulation.

	Target/Background (s-1)
	Contrast=4	Contrast=3	Contrast=2
$K_{p e}$	0.012/0.003	0.009/0.003	0.006/0.003
$K_{e p}$	0.004/0.001	0.003/0.001	0.002/0.001
$K_{p}$	0.025/0.025

The pharmacokinetic-rate images, $K_{p e} (r)$ and $K_{e p} (r)$ , are estimated using the proposed direct strategy by setting the noise covariance matrices in the filter to $Q = 10^{- 7} I$ and $R = 1 \times 10^{- 4} I$ , and the initial pharmacokinetic-rates to 0, and are compared with the results reconstructed by the indirect strategy.¹⁴^,²⁰ To simplify the procedure, an a priori value of $K_{p}$ is empirically determined by fitting the time course of the normalized Born measurement at a fixed detection position, $I_{n b} (k)$ , and a single exponential function $B \exp (- K_{p} Δ T k)$ .³¹ Figure 3 shows the estimated results for three contrasts of $K_{p e}$ and $K_{e p}$ including the images and their X-profiles by the indirect strategy and the direct strategy. The black circles indicate the ideal location and size of the targets. As shown in Fig. 3, a good agreement between the true and the estimated images is observed in terms of the localization and size of the target. It can be found from the X-profiles that the direct strategy outperforms the indirect one in the estimation accuracy. Three metrics are introduced for quantitative assessment of the accuracy and sensitivity performances, referred to as mean square error (MSE), contrast-to-noise ratio (CNR), and quantitativeness ratio (QR).³²

Eq. (6)

{\begin{cases} MSE = \frac{1}{N} \sum_{n = 1}^{N} {[\hat{X} (r_{n}) - X (r_{n})]}^{2} \\ CNR = \frac{\sqrt{\frac{1}{N_{t}} \sum_{n = 1}^{N_{t}} [\hat{X} (r_{n}) - X_{b g}]}}{\sqrt{\frac{1}{N_{b g}} \sum_{n = 1}^{N_{b g}} [\hat{X} (r_{n}) - X_{b g}]}} \\ QR = \frac{\max_{n} [\hat{X} (r_{n})]}{\max_{n} [X (r_{n})]} \end{cases},

where

\hat{X} (r_{n})

and

X (r_{n})

are the estimated and the true pharmacokinetic parameters at the

n

’th node, respectively;

N_{t}

and

N_{b g}

are the node numbers in the target and the background regions, respectively; and

X_{b g}

is the mean value of the true pharmacokinetic parameter in the background area. Table 3 lists the three metrics calculated for the above simulations, from which it is seen that the direct strategy achieves a lower MSE, a higher CNR, and a higher QR than the indirect one.

Fig. 3

Estimated results for three contrasts of $K_{p e}$ (up) and $K_{e p}$ (bottom): (a) images of the pharmacokinetic-rates by the indirect strategy, (b) images of the pharmacokinetic-rates by the direct strategy, and the X-profiles for the scenario of, and (c) $contrast = 2$ , (d) $contrast = 3$ , and (e) $contrast = 4$ of the two strategies. The black circles indicate the ideal location and size of the targets.

Table 3

Metrics of the estimations by the two strategies with three contrasts.

	Contrast	Indirect/Direct
	Contrast	MSE(×10−7)	CNR	QR
$K_{p e}$	2	92.34/25.99	0.99/1.77	0.024/0.77
	3	105.53/27.82	0.98/1.56	0.020/0.67
	4	126.96/26.76	0.99/1.84	0.012/0.56
$K_{e p}$	2	10.62/7.33	0.93/1.40	0.081/0.77
	3	12.03/5.17	0.86/1.66	0.093/0.85
	4	13.49/8.24	0.70/1.26	0.17/0.64

In addition to the pharmacokinetic-rates, the time-courses of the EES and plasma concentration images, $C_{e} (k) = {[C_{e} (r_{1}, k), \dots, C_{e} (r_{N}, k)]}^{T}$ and $C_{p} (k) = {[C_{p} (r_{1}, k), \dots, C_{p} (r_{N}, k)]}^{T}$ , are also estimated in the procedure, giving the total concentration image of $C_{T} (k) = C_{e} (k) + C_{p} (k)$ . Figure 4 shows the estimated average time-courses of ICG compartmental concentrations and their MSEs in the target and background regions for the scenario of $contrast = 3$ . It is expected that the direct estimations of the average compartmental concentrations ${\bar{C}}_{e} (k)$ and ${\overline{C}}_{p} (k)$ , (therefore, the average total concentration ${\overline{C}}_{T} (k)$ , converge faster than the indirect estimations in both the target and background regions. This observation presents a reinforced validation of the superiority of the direct estimation to the indirect one.

Fig. 4

Estimated average time-courses of ICG concentrations with ${\bar{C}}_{e} (k)$ (top), ${\bar{C}}_{p} (k)$ (middle), and ${\overline{C}}_{T} (k)$ (bottom) and their MSEs for the scenario of $contrast = 3$ : (a) in the target and (b) background regions.

3.2.

Phantom Experiments

A pilot phantom experiment is conducted using the lab-equipped CW-DFT system of CT-analogous scanning mode, as shown in Fig. 5. The system uses a 780-nm diode laser (MDL300, PicoQuant), specifically suitable for the ICG excitation. The laser beam is coupled into a source fiber with a core diameter of $300 μ m$ and a numerical aperture (NA) of 0.22, then collimated to impinge the boundary of the phantom. The transmitted light is collected by eight ( $D = 8$ ) detection fibers with a $500 - μ m$ core diameter and $NA = 0.37$ , and coupled into a $8 \times 4$ fiber-optic switch with its output collimated again for normal incidence to four motorized filter wheels (FW102C, Thorlabs) that house a single-bandpass interference filter (ICG-B Emitter, Semrock) each with a center wavelength of the passband 830 nm for the detection of the ICG emission. The filtered beams finally enter into four ( $N_{d} = 4$ ) PMTs photon-counting heads (H8259-02, Hamamatsu, Japan) and are resolved by a four-channel counter at an integration time (i.e., sampling period) of $Δ T = 1 s$ . By rotating the imaging chamber at an angular interval of 22.5 deg, four ( $P = 4$ ) complete sets of the 360 deg projections can be acquired.

Fig. 5

(a) Experimental setup. (b) Phantom and imaging plane.

A liquid phantom is used as shown in Fig. 5(b), which is made from a cylindrical polyformaldehyde chamber of 30-mm diameter filled with 1%-Intralipid as the background and a cylindrical tube of 12-mm diameter filled with a mixture of 1%-Intralipid solution and ICG (GE Healthcare) with an initial concentration of $C_{0} = 8 μ M$ as the target. To simulate the simple one-compartment kinetic process in the target region, 1%-Intralipid fresh solution and the base solution are pumped into and out of the target tube, respectively, by two pumps (BT100-02, China) simultaneously working at the same flow velocity $v_{p}$ , keeping a constant target volume of $V_{0} = 5 mL$ .

The system matrix in Eq. (3) reduces to $K [θ (r)] = τ (r) = \exp (- K_{p} Δ T)$ with $K_{p} = - v_{p} / V_{0}$ , for the one-compartment model in this scenario. Three sets of the flow velocities, i.e., $v_{p} = 0.0053, 0.021$ , and $0.042 mL / s$ , are employed in the validation to simulate the different kinetic processes. Figure 6 shows the $K_{p}$ images and target concentration time-courses estimated by the proposed direct scheme. It can be found that, for all the three cases, the $K_{p}$ is estimated with reasonably small errors in view of the localization and size of the target, as well as approximate proportions to the expected values, as shown in Fig. 6(c). Although the estimations of the average target concentrations deviate from the expected ones, the different kinetic processes are explicitly distinguished.

Fig. 6

Experimental results estimated by the direct strategy: (a) the images of $K_{p}$ , (b) the X-profiles of $K_{p}$ , (c) the QRs of $K_{p}$ , and (d) the average target concentration time-courses.

4. Discussions and Conclusions

In the above validations, we have presented the quantitative comparison between the proposed direct strategy and the conventional indirect one for different kinetic contrasts, based on a cost-effective system configuration using four parallel detection channels. Essentially, optimization of the illumination-detection configuration or the sampling period in the direct strategy within the adaptive-EKF scheme might be feasible for enhancement of the methodology performance, and more importantly, for simplification of the instrumentation complexity, with a criterion of retaining a constant total time of the kinetic process.

4.1.

Illumination-Detection Configuration

The performances of the proposed direct scheme are simulatively investigated using four different illumination-detection configurations, as listed in Table 4. All the configurations work with the aforementioned scenario with kinetic contrast of three, the sampling interval of $Δ T = 2.5 s$ and the SNRs of 50 dB and 40 dB for the excitation and emission, respectively. Figure 7 shows the estimated $K_{p e}$ - and $K_{e p}$ -images, from which it can be found that, although the three assessment metrics calculated in Table 4 quantitatively present small differences among the four reconstructions with Case 1 being overall optimized and Case 3 minimally quantitative due to complex reasons that are worthy of future identification, the images of all cases are qualitatively of the same quality. Nevertheless, it is favorable to conclude that, without evident sacrifice of the performances, the configuration in Case 2 can implement a significant reduction in the instrumentation complexity and expense by omitting the multiplexing function.

Fig. 7

Estimated results by the direct strategy with different illumination-detection configurations in Table 4: (a) images of $K_{p e}$ (top) and $K_{e p}$ (bottom), (b) X-profiles of $K_{p e}$ , and (c) X-profiles of $K_{e p}$ .

Table 4

Metrics of the estimations by the direct strategy with different illumination-detection configurations.

Case	S	D	Nd	P	Kpe/Kep
Case	S	D	Nd	P	MSE(×10−7)	CNR	QR
1	16	8	8	16	27.02/4.83	1.64/2.16	0.69/0.95
2	32	4	4	8	29.09/4.97	1.49/2.04	0.68/0.97
3	32	8	8	8	28.23/4.89	1.29/1.83	0.59/0.85
4	8	8	8	32	29.99/5.31	1.40/1.76	0.65/0.90

4.2.

Sampling Period

Physically, given the efficiency of the photoelectric conversion, the SNR of the photon counting signal is proportional to the square root of the integration time. This integration time might be properly chosen for a balance between the noise level and temporal resolution of the sampling, both observably influencing the performances. Assuming fixed SNRs of $α_{x} = 50 dB$ and $α_{m} = 40 dB$ for the reference sampling period of 2.5 s, the SNRs for an arbitrary period $Δ T$ can be simply obtained: $α_{x / m} (Δ T) = α_{x / m} (2.5 s) + 20 \log_{10} (\sqrt{Δ T / 2.5 s})$ . Three sampling periods, as listed in Table 5, are employed for the investigations, which are conducted for the scenario with kinetic contrast of three and the illumination-detection configuration used as Case 1 aforementioned. The pharmacokinetic-rate images, $K_{p e}$ and $K_{e p}$ , estimated by the direct strategy, for the three periods ( $Δ T = 5, 2.5$ , and 1.25 s), are shown in Fig. 8, and the assessment metrics are compared in Table 5. It is found that, although the smaller $Δ T$ generates a larger dataset, the resultant MSE for $Δ T = 2.5 s$ is smallest in the three cases. This manifests a sampling period of 2.5 s as the optimal one under the assumed noise levels.

Fig. 8

Estimated results by the direct strategy with different sampling periods in Table 5: (a) images of $K_{p e}$ (top) and $K_{e p}$ (bottom), (b) X-profiles of $K_{p e}$ , and (c) X-profiles of $K_{e p}$ .

Table 5

Metrics of the estimations by the direct strategy with different sampling period.

ΔT(s)	αx/αm (dB)	P	Kpe/Kep
ΔT(s)	αx/αm (dB)	P	MSE(×10−7)	CNR	QR
5	53/43	8	39.43/3.80	1.04/1.95	0.70/0.89
2.5	50/40	16	27.02/4.83	1.64/2.16	0.69/0.95
1.25	47/37	32	55.62/4.95	2.13/2.72	0.99/1.22

In conclusion, both the simulations and the phantom experiments reveal that the proposed direct strategy estimates the compartmental concentrations and the pharmacokinetic-rates with a fair quantitative and localization accuracy. This strategy relaxes the restrictions on the data-acquisition mode and the system configuration, achieving a highly-sensitive and cost-effective detection based on the photon-counting DFT with CT-scanning mode. Future work will focus on the in vivo experiments of small animals and exploring the optimization in term of the illumination-detection configuration and the sampling period, and so on.

Acknowledgments

The authors acknowledge the funding supports from the National Natural Science Foundation of China (81271618, 81371602, 61475115, 61475116, 81401453, 61575140, and 81571723), and Tianjin Municipal Government of China (13JCZDJC28000, 14JCQNJC14400, and 15JCZDJC31800).

References

1.

K. M Tichauer et al., “Quantitative in vivo cell-surface receptor imaging in oncology: kinetic modeling and paired-agent principles from nuclear medicine and optical imaging,” Phys. Med. Biol., 60 R239 –R269 (2015). http://dx.doi.org/10.1088/0031-9155/60/14/R239 Google Scholar

2.

H. Maeda, “The enhanced permeability and retention (EPR) effect in tumor vasculature: the key role of tumor-selective macromolecular drug targeting,” Adv. Enzyme Regul., 41 189 –207 (2001). http://dx.doi.org/10.1016/S0065-2571(00)00013-3 AEZRA2 0065-2571 Google Scholar

3.

N. Kosaka et al., “Near infrared fluorescence-guided real-time endoscopic detection of peritoneal ovarian cancer boules using intravenously injected indocyanine green,” Int. J. Cancer, 129 1671 –1677 (2011). http://dx.doi.org/10.1002/ijc.v129.7 Google Scholar

4.

L. A. Bauer, Applied Clinical Pharmacokinetics, McGraw-Hill, New York (2008). Google Scholar

5.

X. Intes et al., “In vivo continuous-wave optical breast imaging enhanced with indocyanine green,” Med. Phys., 30 1039 –1047 (2003). http://dx.doi.org/10.1118/1.1573791 MPHYA6 0094-2405 Google Scholar

6.

X. Liu et al., “Imaging of indocyanine green perfusion in mouse liver with fluorescence diffuse optical tomography,” IEEE Trans. Biomed. Eng., 58 2139 –2143 (2011). http://dx.doi.org/10.1109/TBME.2011.2135858 IEBEAX 0018-9294 Google Scholar

7.

B. Alacam et al., “Extended Kalman filtering for the modeling and analysis of ICG pharmacokinetics in cancerous tumors using NIR optical methods,” IEEE Trans. Biomed. Eng., 53 1861 –1871 (2006). http://dx.doi.org/10.1109/TBME.2006.881796 IEBEAX 0018-9294 Google Scholar

8.

B. Alacam et al., “Pharmacokinetic-rate images of indocyaninegreen for breast tumors using near-infrared optical methods,” Phys. Med. Biol., 53 837 –859 (2008). http://dx.doi.org/10.1088/0031-9155/53/4/002 PHMBA7 0031-9155 Google Scholar

9.

D. J. Cuccia et al., “In vivo quantification of optical contrast agent dynamics in rat tumors by use of diffuse optical spectroscopy with magnetic resonance imaging coregistration,” Appl. Opt., 42 2940 –2950 (2003). http://dx.doi.org/10.1364/AO.42.002940 APOPAI 0003-6935 Google Scholar

10.

A. B. Milstein, K. J. Webb and C. A. Bouman, “Estimation of kinetic model parameters in fluorescence optical diffusion tomography,” J. Opt. Soc. Am. A, 22 1357 –1368 (2005). http://dx.doi.org/10.1364/JOSAA.22.001357 JOAOD6 0740-3232 Google Scholar

11.

G. Zhang et al., “A direct method with structural priors for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography,” IEEE Trans. Biomed. Eng., 61 986 –990 (2014). http://dx.doi.org/10.1109/TBME.2013.2292714 IEBEAX 0018-9294 Google Scholar

12.

G. Zhang et al., “Full-direct method for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography,” Appl. Phys. Lett., 106 081110 (2015). http://dx.doi.org/10.1063/1.4913690 APPLAB 0003-6951 Google Scholar

13.

G. Zhang et al., “Bayesian framework based direct reconstruction of fluorescence parametric images,” IEEE Trans. Med. Imaging, 34 1378 –1391 (2015). http://dx.doi.org/10.1109/TMI.2015.2394476 ITMID4 0278-0062 Google Scholar

14.

B. Alacam and B. Yazici, “Direct reconstruction of pharmacokinetic-rate images of optical fluorophores from NIR measurement,” IEEE Trans. Med. Imaging, 28 1337 –1353 (2009). http://dx.doi.org/10.1109/TMI.2009.2015294 ITMID4 0278-0062 Google Scholar

15.

V. Ntziachristos et al., “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med., 8 757 –761 (2002). http://dx.doi.org/10.1038/nm729 1078-8956 Google Scholar

16.

A. B. Milstein et al., “Fluorescence optical diffusion tomography,” Appl. Opt., 42 3081 –3094 (2003). http://dx.doi.org/10.1364/AO.42.003081 APOPAI 0003-6935 Google Scholar

17.

J. Baeten et al., “Development of fluorescent materials for diffuse fluorescence tomography standards and phantoms,” Opt. Express, 15 8681 –8694 (2007). http://dx.doi.org/10.1364/OE.15.008681 OPEXFF 1094-4087 Google Scholar

18.

L. Zhang et al., “Three-dimensional strategy for time-domain fluorescence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express, 16 7214 –7223 (2008). http://dx.doi.org/10.1364/OE.16.007214 OPEXFF 1094-4087 Google Scholar

19.

F. Gao et al., “Simultaneous fluorescence yield and lifetime tomography from time-resolved transmittances of small-animal-sized phantom,” Appl. Opt., 49 3163 –3172 (2010). http://dx.doi.org/10.1364/AO.49.003163 APOPAI 0003-6935 Google Scholar

20.

X. Wang et al., “Performance enhancement of pharmacokinetic diffuse fluorescence tomography by use of adaptive extended Kalman filtering,” Comput. Math. Methods Med., 2015 1 –15 (2015). http://dx.doi.org/10.1155/2015/739459 Google Scholar

21.

X. Liu et al., “4-D reconstruction for dynamic fluorescence diffuse optical tomography,” IEEE Trans. Med. Imaging, 31 2120 –2132 (2012). http://dx.doi.org/10.1109/TMI.2012.2213828 ITMID4 0278-0062 Google Scholar

22.

X. Zhang et al., “Reconstruction of fluorophore concentration variation in dynamic fluorescence molecular tomography,” IEEE Trans. Biomed. Eng., 62 138 –144 (2015). http://dx.doi.org/10.1109/TBME.2014.2342293 IEBEAX 0018-9294 Google Scholar

23.

X. Wang et al., “A high-sensitive diffuse fluorescence tomography system with CT-analogous scanning mode,” Proc. SPIE, 82160V (2012). http://dx.doi.org/10.1117/12.905355 Google Scholar

24.

V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett., 26 893 –895 (2001). http://dx.doi.org/10.1364/OL.26.000893 Google Scholar

25.

S. Sourborn, “A tracer-kinetic field theory for medical imaging,” IEEE Trans. Med. Imaging, 33 935 –946 (2014). http://dx.doi.org/10.1109/TMI.2014.2300450 ITMID4 0278-0062 Google Scholar

26.

S. P. Sourbron and D. L. Buckley, “Tracer kinetic modeling in MRI: estimating perfusion and capillary permeability,” Phys. Med. Biol., 57 R1 –R33 (2012). Google Scholar

27.

D. H. Anderson, “Compartmental modeling and tracer kinetics,” Lecture Notes in Biomathematics, 50 Springer-Verlag, Berlin, Germany (1983). Google Scholar

28.

J. A. Jacquez, Compartmental Analysis in Biology and Medicine: Kinetics of Distribution of Tracer-Labeled Materials, Elsevier, New York (1972). Google Scholar

29.

M. Boutayeb and D. Aubry, “A strong tracking extended Kalman observer for nonlinear discrete-time systems,” IEEE Trans. Autom. Control, 44 1550 –1556 (1999). http://dx.doi.org/10.1109/9.780419 IETAA9 0018-9286 Google Scholar

30.

L. Ozbek and M. Efe, “An adaptive extended Kalman Filter with application to compartment models,” Commun Stat. Simul. Comput., 33 145 –158 (2004). http://dx.doi.org/10.1081/SAC-120028438 CSSCDB 0361-0918 Google Scholar

31.

T. Desmettre, J. M. Devoisselle and S. Mordon, “Fluorescence properties and metabolic features of indocyanine green (ICG) as related to angiography,” Surv. Ophthalmol., 45 15 –27 (2000). http://dx.doi.org/10.1016/S0039-6257(00)00123-5 SUOPAD 0039-6257 Google Scholar

32.

L. Wu et al., “Enhancement of fluorescence molecular tomography with structural-prior-based diffuse optical tomography: combating optical background uncertainty,” Appl. Opt., 53 6970 –6982 (2014). http://dx.doi.org/10.1364/AO.53.006970 APOPAI 0003-6935 Google Scholar

Biography

Xin Wang is a PhD candidate at Tianjin University. She received her BS and MS degrees in biomedical engineering from Tianjin University in 2010 and 2013, respectively. Her current research interests include fluorescence tomography and pharmacokinetic imaging, and she has authored more than five journal papers.

Biographies for the other authors are not available.

Citation Download Citation

Xin Wang, Yanqi Zhang, Limin Zhang, Jiao Li, Zhongxing Zhou, Huijuan Zhao, and Feng Gao "Direct reconstruction in CT-analogous pharmacokinetic diffuse fluorescence tomography: two-dimensional simulative and experimental validations," Journal of Biomedical Optics 21(4), 046007 (19 April 2016). https://doi.org/10.1117/1.JBO.21.4.046007

Published: 19 April 2016

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 10 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Fluorescence tomography

Data modeling

Tissues

Signal to noise ratio

Instrument modeling

Performance modeling

Digital filtering

1.

Introduction

2.

Methods

Eq. (1)

Fig. 1

2.1.

Kinetics-to-Measurement Mapping

Eq. (2)

Eq. (3)

Eq. (4)

2.2.

Adaptive-Extended Kalman Filtering

Table 1

3.

Validations

3.1.

Simulative Investigations

Eq. (5)

Fig. 2

Table 2

Eq. (6)

Fig. 3

Table 3

Fig. 4

3.2.

Phantom Experiments

Fig. 5

Fig. 6

4.

Discussions and Conclusions

4.1.

Illumination-Detection Configuration

Fig. 7

Table 4

4.2.

Sampling Period

Fig. 8

Table 5

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years