Abstract
Due to both the physiological and morphological differences in the vascularization between healthy and diseased tissues, pharmacokinetic diffuse fluorescence tomography (DFT) can provide contrast-enhanced and comprehensive information for tumor diagnosis and staging. In this regime, the extended Kalman filtering (EKF) based method shows numerous advantages including accurate modeling, online estimation of multiparameters, and universal applicability to any optical fluorophore. Nevertheless the performance of the conventional EKF highly hinges on the exact and inaccessible prior knowledge about the initial values. To address the above issues, an adaptive-EKF scheme is proposed based on a two-compartmental model for the enhancement, which utilizes a variable forgetting-factor to compensate the inaccuracy of the initial states and emphasize the effect of the current data. It is demonstrated using two-dimensional simulative investigations on a circular domain that the proposed adaptive-EKF can obtain preferable estimation of the pharmacokinetic-rates to the conventional-EKF and the enhanced-EKF in terms of quantitativeness, noise robustness, and initialization independence. Further three-dimensional numerical experiments on a digital mouse model validate the efficacy of the method as applied in realistic biological systems.
1. Introduction
In diffuse fluorescence tomography (DFT) regime, pharmacokinetic imaging means a dynamic modality that ultimately acquires the spatially varying pharmacokinetic parameters of administrated fluorescent agent in tissue [1]. Among all the commercially available fluorescent agents, only indocyanine green (ICG) is approved for human use by the U.S. Food and Drug Administration. As a blood pooling agent, ICG has evidently distinct delivery behavior between cancerous tissue and normal tissue due to the proliferation of the “leaky” angiogenetic microvessels [2]. Therefore, the pharmacokinetic-DFT of ICG potentially provides contrast- and specificity-enhanced information for tumor diagnosis, malignancy staging, treatment monitoring, and drug-delivery assessment, as compared to the static modality that only discloses the temporally averaged fluorophore concentration image [3].
The current image reconstruction methods for the pharmacokinetic-DFT can be categorized into explicit and implicit schemes, both combining the DFT principle with pharmacokinetic analysis. In the explicit scheme the pharmacokinetic images are calculated in a voxel-by-voxel fashion by fitting a kinetic model to the reconstructed temporal sequence of the fluorophore concentration [4], while in the implicit scheme the pharmacokinetic images are directly reconstructed as a whole by incorporating a deterministic kinetic model to the inversion procedure [5] or by expressing the kinetics-to-measurement map in the extended Kalman filtering (EKF) procedure [6]. Although it is demonstrated that the implicit scheme substantially improves accuracy and robustness in pharmacokinetic parameter estimation, the methodology is theoretically complex and computationally costly due to the further degradation of ill-posedness and increased unknowns in the inversion, thus limiting its application to two-dimensional (2D) scenarios [5, 6]. In this sense, the explicit scheme is universally applicable and relatively easy-to-implement but needs to be enhanced for its performance and robustness in both the concentration reconstruction and the parameter estimation.
Compartmental modeling is a well-established approach to the pharmacokinetic analysis. This method describes the concentration dynamics as a result of fluorophore exchange among kinetically distinct compartments using a set of coupled ordinary differential equations (ODEs) with its coefficients representing the exchange rates, referred to as the pharmacokinetic-rates. A biexponential-curve-fitting method based on the two-compartment model has been proposed to demonstrate the feasibility of the ICG pharmacokinetics in tumor diagnosis [4, 7–9]. Alacam and Yazici have addressed an EKF framework for estimating the pharmacokinetic-rates of highly nonlinear nature and validated the sufficiency of the two-compartment model in describing ICG kinetics [3]. Furthermore, an EKF study using the two-compartment model has indicated the superiority of the pharmacokinetic-rate images to the bulk rates of entire breast in cancer diagnosis [10]. In comparison with the curve-fitting techniques, the EKF-based method has a number of advantages: (1) the EKF regards the ICG dynamics as an evolutionary stochastic process and therefore provides a better fit to the fluorophore metabolism than the exponential-curve-based model; (2) the EKF can readily accommodate spatiotemporal priors of the kinetic parameters; (3) the EKF can achieve real-time estimation of a complete set of the pharmacokinetic parameters including pharmacokinetic-rates and concentrations in different compartments; and (4) the method is universally applicable to any optical fluorophore [3, 6].
However, the performance of the conventional-EKF highly hinges on the exact prior knowledge about the initial states, that is, the expectations and covariances of the compartmental concentrations and pharmacokinetic parameters, which are always inaccessible in practice. The inappropriate initialization can cause the degeneration even the divergence of the EKF [11–13]. Ozbek et al. have presented an enhanced-EKF scheme to improve the accuracy of the EKF analysis [14, 15], where a forgetting-factor is introduced to compensate the inaccuracy of the initial states and to emphasize the effect of the current data. This strategy in principle better stabilizes the online estimation compared with the conventional-EKF [16]. Nevertheless, a constant forgetting-factor is intrinsically suboptimal, for example, the larger forgetting-factor results in a larger Kalman gain means that the filter can quickly adapt to a new situation but also that it is sensitive to random errors. By nature, adopting a variable forgetting-factor that adaptively balances between the convergence and robustness will further improve the performance of the EKF.
We herein propose an adaptive-EKF for DFT-based pharmacokinetic imaging, where the forgetting-factor is updated at each recursive stage, by contrasting the calculated innovation covariance and the estimated one. For the methodology to be universally applicable, the algorithm is implemented within the framework of the explicit scheme that firstly reconstructs the time-course of ICG concentration with the conventional-DFT and then accordingly estimates the pharmacokinetic-rate images using the adaptive-EKF approach for the two-compartment model. Simulation results of a two-dimensional circular model suggest that the proposed adaptive-EKF can obtain preferable pharmacokinetic-rate images to both the conventional- and enhanced-EKF in terms of quantitativeness, noise robustness, and initialization independence. Further numerical experiments on a three-dimensional (3D) digital mouse model validate the feasibility and efficacy of the method as applied in 3D complex biological systems.
2. Theory
2.1. Two-Compartment Model
Compartmental modeling assumes that a biological system is conceptually (not geometrically) divided into a series of compartments, each representing a well-mixed space of similar tissues within which the fluorophore is uniformly distributed and its concentration changes as a result of the agent exchange among the compartments. Mathematically, the kinetic changes of the compartmental concentrations are governed by a collection of coupled ODEs, ultimately resting on the principle of mass conservation [17–20]. In the two-compartment model, tissue is composed of plasma and the extracellular-extravascular space (EES), as shown in Figure 1, where and denote the tissue concentrations of ICG in plasma and EES, that is, the numbers of ICG molecules in the two compartments relative to the total tissue volume, respectively; , , and are pharmacokinetic-rates describing the ICG leakage into and the drainage out of the EES, as well as the ICG elimination from the body through circulatory system, respectively; is the arterial input function (AIF), and is volume flow [17]. With the above notations, the two-compartment model of ICG is described by the following ODEs set:
It is important to distinguish the tissue compartmental concentrations, (), used in the above ODE from the local compartmental concentrations, , with and being the fractions of plasma and EES volumes, respectively [17, 18]. The latter is essentially the ratio of the ICG molecule number in plasma or EES compartment to its fractional volume and has been inappropriately used in the previous works [3, 6, 10, 21].
With the explicit scheme of dynamic DFT, the time course of the total ICG concentration in tissue, , is tomographically reconstructed at discrete time instances, , , where is the sampling period. Here we only consider the permeability of ICG in the EKF process after agent administration. In order to achieve the joint estimation of the pharmacokinetic-rates and the ICG tissue concentrations within the EKF framework, a dynamic model of the parameter vector is additionally appended to the two-compartment model to construct a discrete nonlinear state-space model as follows [3, 6, 17–20]: where is the vector representing the compartmental concentrations with and ; ; , , and are independent zero-mean Gaussian white noise processes with the covariance matrices and and the variance , referred to as the state driving noise vector, parameter driving noise vector, and observation noise, respectively; is the parameter estimation at time ; is the pharmacokinetic-rates-related system matrix for the discrete time two-compartment model as follows:
2.2. The General Framework of EKF
Let and be the estimation of the compartmental concentrations and the parameters at the th step, respectively, let and be the one-step ahead prediction of and , respectively. The general framework of an EKF for the two-compartment model is summarized as a recursive procedure of “Prediction-Gain-Update” [4, 6, 13–16].
(1) Prediction. Consider
(2) Kalman Gain. Considerwhere , is the Jacobian matrix of the nonlinear EKF system, is referred to as the calculated covariance since it equals the covariance of the innovation sequence as the model is accurate [13], is an identity matrix, and is a forgetting-factor. To obtain the Kalman gain, , the one-step ahead prediction of the error-covariance matrix, , is firstly calculated in terms of the estimated error-covariance matrix at the previous step, , and then is updated to reach its estimation at the current step, , for the next recursion. According to the forgetting-factor , the EKF can be classified as the conventional () [4, 6, 10] or the enhanced ( is a constant larger than 1) [14, 15].
(3) Update. Consider
Prior to the recursive process, the EKF can be initialized theoretically for the state, that is, expectations of ICG concentrations and and the error covariance matrix, that is, . In practice, is always experientially chosen as the physiologically relevant values. can be calculated from the initial quantities of intravenous injection. It is emphasized that the filter may degenerate or even diverge with the improper initialization [11–13].
2.3. Adaptive-EKF with a Variable Forgetting-Factor
Theoretically, a linear filter can be considered an optimal one as its innovation sequence is white. According to this criterion, however, the innovation sequence of the EKF is not white due to the linearization error. Nevertheless, the performance of the EKF might still be improved on a condition of temporal independence of the innovation; that is, the autocovariance of the innovation sequence is zero. To approach such a condition, a temporally varying forgetting-factor is introduced as [16]where is referred to as the estimated innovation covariance:The essence of the above equation is to estimate the real innovation covariance by averaging inside a moving estimation window of size [22]. The windowing strategy in the calculation helps balance between the confidences in the “old” observation data and the current data. No optimal criterion is thus far found for selection of the window width and some care must be taken with this issue. In principle, the filtering might fail to mitigate the negative effects of the “old” observation data with a large and, on the contrary, would be adversely affected by sudden noise in the new data as is too small.
It is difficult in practice to acquire the exact nonlinear stochastic equations of the pharmacokinetic system, due to the unavailability of the noise characteristics. With inaccurate initialization of the noise characteristics, the calculated innovation covariance may be lower than the estimated one , leading to a forgetting-factor greater than 1, that is, . The predicted error covariance is then amplified by the forgetting-factor to ameliorate the inexact modeling and emphasize the role of the current data. On the other hand, it must be ensured that to stabilize the filtering process. Since a sharp variation of the ratio of the estimated innovation covariance to the calculated one may occur during the rapidly varying early phase of the ICG kinetic process after injection, due to the inappropriate initialization of the filter or some unaccounted perturbation [2], a logarithmic ratio value of to is used to avoid the rapid local convergence under a large forgetting-factor.
3. Simulative Investigations
Using the pharmacokinetic-rates listed in Table 1, the time-course of the ICG concentration, , is calculated based on the two-compartment kinetic model, that is, (2), that are driven by the zero-mean Gaussian noises, , with signal-to-noise ratio (SNR), , respectively. The pharmacokinetic-rates are assumed to be nearly constant at the clinical-pathologic stage as the measurement is performed, a high SNR of (0.1%) is used to drive the parameter vector in the simulations. In the calculation, we set the initial concentration of ICG to and assume that all the ICG is in the plasma at = 0 s.
With the time-course of the ICG concentration, , we can accordingly calculate the time-varying fluorephore absorption coefficient, , in terms of a linear relationship of , where is the extinction coefficient of ICG, and finally obtain the boundary flux by solving the coupled diffusion equations using the finite element method, as the simulated data for the reconstruction [23–27]. The investigations herein rely on combination of the multichannel photon-counting DFT system custom-made in our lab [21] and the widely-adopted normalized Born formulation for DFT reconstruction [23, 27]. This means that both the excitation and emission SNRs, that is, and , can achieve reasonably high levels of above 30 dB and 20 dB, respectively, and an optically homogeneous background can be used in the forward calculation.
To further obtain the estimated pharmacokinetic parameters and , the DFT-reconstructed ICG concentration, , is analyzed by conventional-, enhanced-, and adaptive-EKF procedures, respectively. The constant forgetting-factor of the enhanced-EKF is set to be 1.1 in the simulations according to [15], while the adaptive-EKF adopts a 7-length () window. The noise covariance matrices of the state and parameter driving noise vectors in the filter are set to and , and the variance of the observation noise is dependent on the additive noise levels in the measurements [6].
3.1. 2D Circular Phantom
A circular phantom with 30 mm diameter is used that embeds a circular tumor-emulating region of 6 mm diameter, as shown in Figure 2. Its absorption coefficient and reduced scattering coefficient are set to and for both the excitation and emission wavelengths. The domain is discredited by a finite element mesh with 721 nodes and 1350 elements (triangles) for both the forward and inverse calculations and sampled by 16 coaxial source-detector optodes that are placed around the phantom with equal spacing. To acquire a complete dataset for the 16 × 16 source-detector combinations, the 16 detectors collect the photons in parallel as 16 sources illuminate the surface in serial. The measurement is repeatedly conducted for 720 s (12 minutes) at a sampling period of = 10 s, generating 72 datasets of the time-course for the pharmacokinetic estimation.
For the pharmacokinetic-rates, and are set to represent three target-to-background contrasts in agreement with the different tumor pathological staging, while is assumed to be homogeneous and known throughout the simulations, as shown in Table 1.
To demonstrate the above dynamic DFT procedure, Figure 3 illustrates the time-course of the model-simulated and DFT-reconstructed ICG average concentrations in the background and target regions, and , as well as the time-course of the true and reconstructed contrasts, for the scenario of contrast = 2 in Table 1. Here we define the contrast as the ratio of the average target variation to the average background [28]. For the reference, the interim DFT-reconstructed yield-images and their X-profiles are also given at six time instants of 100 s to 600 s, which are generated from the simulated time-course of ICG concentration with , , , and . It is found that the reconstructed concentration features an increased underestimation with increased true contrast during the kinetic process: an adversity originates from the underestimation of high yield contrasts in DFT-reconstruction.
(a)
(b)
(c)
3.1.1. Quantitativeness
The images of the pharmacokinetic-rates, and , estimated by the conventional-, enhanced-, and adaptive-EKFs, for the three sets of the true target-to-background contrasts in Table 1, are illustrated in Figures 4(a) and 4(b), respectively. In the filtering process, the variance of the observation noise is chosen as , and the initial pharmacokinetic-rates are set to be same for the three filters, that is, and .
(a)
(b)
(c)
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 in terms of the X-profiles that the adaptive-EKF outperforms the conventional-EKF and enhanced-EKF in estimation accuracy. For quantitative assessment of the method, the quantitativeness ratio (QR), defined as the ratio of the estimated contrast to the true one of a pharmacokinetic-rate, as shown in Figure 4(c). The much higher QRs of the - and -images achieved by the adaptive-EKF exhibit the enhanced ability of the method to dynamically compensate the initial inaccuracies. This feature is further quantified with the time-course of the mean square error (MSE) between the estimated and true compartmental concentration images, defined aswhere ; is the number of the nodes in the region of interest; and are the estimated and the true compartmental concentrations at the th node, respectively. Figure 5 shows the estimated average time-course of ICG compartmental concentrations and their MSEs in the target and background areas. It is observed, for the three pharmacokinetic contrasts, that the average in the target area estimated by adaptive-EKF converges faster to the true value and realizes a smaller MSE than the conventional and the enhanced ones. For estimation, the three filters exhibit almost the same performance according to the figure.
(a)
(b)
(c)
3.1.2. Noise Robustness
To evaluate the noise robustness of the three filters, the estimations are conducted for the pharmacokinetic contrast of 2, with different levels of the state driving noise and measurement noise, as shown in Table 2. The initial pharmacokinetic-rates are set to be same for the three filters, that is, and .
Firstly, to compare the robustness of the three filters to the state driving noise, the data is generated with , 30 dB, and 40 dB, by fixing reasonably low levels of the measurement noise: and , that is, Case 1, Case 2, and Case 3 in Table 2. A fixed observation noise variance of is used in the filtering process. Figure 6 contrasts the resultant QRs of and for the three cases. It is seen that the QRs achieved by the adaptive-EKF are much higher than the other two schemes.
(a)
(b)
Next, the robustness of the filters to the measurement noise is investigated with a fixed high of 40 dB but varying levels of the measurement noise, that is, Cases 4–6 in Table 2. This time the variance of the observation noise, , in the filtering process is set to , , , and for Cases 3–6, respectively. The estimated - and -QRs are shown in Figure 7, from which it is again found that the much higher QRs are achieved by the adaptive-EKF as compared to the other two schemes.
(a)
(b)
3.1.3. Independence to Initialization
The independence of the initialization of and is quantified by both the QR and the MSE for the pharmacokinetic contrast of 2, with the same noise levels as in Case 3. The MSE measure is introduced to assess the bias of the estimation and calculated in the whole domain. The variance of the observation noise in the filters is set to for all the cases. The initial pharmacokinetic-rates, and , are identically set for the three filters, in proportion to the true background ones with a varying factor from 0 to 2. The two measures are illustrated in Figure 8, from which it is clear that the adaptive-EKF achieves significantly larger QR and smaller MSE than the others regardless of the deviations of the initial and from their true backgrounds.
(a)
(b)
3.2. 3D Digital Mouse Model
To evaluate the performance of the proposed adaptive-EKF in small animal scenarios, simulated data is generated using the two-compartment model, on a 3D digital mouse atlas (Digimouse) [29]. To facilitate the forward calculation, the mouse model is assumed to be embedded into a cylindrical chamber of 15 mm radius and a 35 mm height filled with the matching fluid, as shown in Figure 9. The cylindrical domain is discretized into 25956 nodes and 47250 elements (prisms) for use with the finite element method. A cylindrical tumor target of 2.5 mm radius and of 6 mm length is placed in the liver with its center at mm, mm, and mm, as shown in Figure 9(b). Five imaging planes along the height (-axis) at mm, −4 mm, 0 mm, 4 mm, and 8 mm are arranged for data acquisition, with each installing 32 coaxial source-detector optodes around the phantom at equal spacing. The optical parameters of the organs in the digital mouse are listed in Table 3, including the absorption and the reduced scattering coefficients, and , at both the excitation () and emission () wavelengths [30]. The background and target pharmacokinetic-rates are set to those for the case of contrast = 2 in Table 1. The levels of the various noises are set according to , , , and . In the filtering process, the variance of the observation noise is chosen as , and the initial pharmacokinetic-rates are identically set for the three filters, that is, and . The top- and side-view images of the estimated pharmacokinetic-rates are shown at mm and mm in Figure 10, respectively. Analogous to the 2D scenarios, the proposed adaptive-EKF greatly improves the quantitativeness of the pharmacokinetic estimations as compared to the conventional- and enhanced-EKF, in terms of the X-profiles ( mm and mm). The results indicate the prospects of the proposed method in ICG kinetics study of diseased mouse models in vivo.
(a)
(b)
(a)
(b)
(c)
4. Discussions and Conclusions
It is clearly seen from the estimated pharmacokinetic-rate images that the target sizes are slightly overestimated and biased, even with the adaptive-EKF. These two defects and also the quantitativeness remain to be further improved. In practice, a successful imaging of the pharmacokinetic-rates by the explicit methods, such as the EKFs described here, is dependent on two crucial factors: one is the fidelity of the DFT reconstruction, that is, how to reconstruct the ICG concentration that approaches the realistic one and another is the effectiveness of the EKF process, that is, how to accurately extract the pharmacokinetic parameters, and , from the DFT-reconstructed concentration .
The former is restrained by the severe ill-posedness of the DFT inversion, normally resulting in a decreasing quantitativeness with the increasing contrast [25–27]. The distorted time-course of ICG concentration inevitably leads to a deviation in the estimation of the pharmacokinetic-rates. To demonstrate the first effect, Figure 11 compares the estimated images and X-profiles of the pharmacokinetic-rates by the adaptive-EKF from the model-simulated and DFT-reconstructed ICG concentrations, that is, and , respectively, for the case of contrast = 2 in Table 1 with a SNR setting of , , , and . It is obvious that the accuracy of the estimation is significantly improved with the model-simulated concentration by bypassing the DFT process. Nevertheless, the DFT process is indispensable due to the inaccessibility of the ICG concentration in tissue. Therefore, the introduction of more advanced DFT reconstruction methodologies, for example, a prior knowledge guided scheme or an efficient explicit implementation, is necessarily requested for enhancement.
The latter significantly depends on the selection of the noise covariances, that is, , , and , in the EKF process. In practice, these noise covariances vary and are unavailable during the measurement period and are suboptimally assumed to be constant in this work. Although the adaptive-EKF can compensate this inaccurate assumption, the prior knowledge or an adaptive updating of the noise covariance may further improve the performance of the estimation.
In conclusion, an adaptive-EKF is developed based on the two-compartment model, for the enhanced estimation of the pharmacokinetic-rates from the dynamic DFT reconstruction. With introduction of a variable forgetting-factor, the propose scheme can effectively compensate the uncertainties of the initial states and the noise covariances. The simulation results suggest that the adaptive-EKF can obtain preferable pharmacokinetic-rate images than the conventional-EKF and the enhanced-EKF with an improved quantitiveness, noise robustness, and initialization independence. The mouse experiments in vivo are necessary to study the real pharmacokinetic process in our future work.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge the funding supports from the National Natural Science Foundation of China (81271618, 81371602, 61475115, 61475116, and 81401453), the Research Fund for the Doctoral Program of Higher Education of China (20120032110056), and the Tianjin Municipal Government of China (12JCQNJC09400, 13JCZDJC28000, 14JCQNJC14400).