Parameterized level-set based pharmacokinetic fluorescence optical tomography using the regularized Gauss–Newton filter

Omprakash Gottam; Naren Naik; Sanjay Gambhir

doi:10.1117/1.JBO.24.3.031010

10 October 2018 Parameterized level-set based pharmacokinetic fluorescence optical tomography using the regularized Gauss–Newton filter

Omprakash Gottam, Naren Naik, Sanjay Gambhir

Author Affiliations +

Journal of Biomedical Optics, Vol. 24, Issue 3, 031010 (October 2018). https://doi.org/10.1117/1.JBO.24.3.031010

Abstract

Pharmacokinetic tomography is emerging as an important methodology for detecting abnormalities in tissue based upon spatially varying estimation of the pharmacokinetic rates governing the leakage of an injected fluorophore between blood plasma and tissue. We present a shape-based reconstruction framework of a compartment-model based formulation of this dynamic fluorescent optical tomography problem to solve for the pharmacokinetic rates and concentrations of the fluorophore from time-varying log intensity measurements of the optical signal. The compartment-model based state variable model is set up in a radial basis function parameterized level set setting. The state (concentrations) and (pharmacokinetic) parameter estimation problem is solved with an iteratively regularized Gauss–Newton filter in a trust-region framework. Reconstructions obtained using this scheme for noisy data obtained from cancer mimicking numerical phantoms of near/sub-cm sizes show a good localization of the affected regions and reasonable estimates of the pharmacokinetic rates and concentration curves.

1. Introduction

The detection of early cancer requires the capturing of physiological changes that occur in the affected cells before their structure/morphology changes. Tomographic modalities that capture the tell-tale biochemical/physiological signatures of early cancer are nuclear medicine-based schemes, such as positron emission tomography, or optical tomography based ones, such as fluorescence optical tomography (FOT). Two kinds of information that can be inferred about the region of interest from these modalities are the absorption (concentration) distribution of introduced markers (that attach themselves to the affected regions), and the rates at which the markers enter and leave the regions of interest.¹^–¹⁴

Pharmacokinetic rates are the parameters that govern the passage of the markers across notional compartments in the body, such as blood-plasma and tissue ones.¹⁵^–¹⁷ These rates have been found to bear the signatures of abnormalities in the tissue,¹^,⁴^,⁹^,¹⁰^,¹⁸^–²² thus yielding a potentially powerful tool of early diagnosis. In FOT, pharmacokinetic-rate reconstructions have been proposed in a compartment model setting, to detect oncological abnormalities in tissue,⁹^,¹⁰^,¹⁹^–²² where a state variable model is constructed with fluorophore-concentrations being the states and the pharmacokinetic rates and volume fractions being the parameters. In Ref. 21, the concentrations at all points in the domain are reconstructed at each time instant from a linearized (Rytov) approximation-based inversion of the complex log-intensity data at each instant. Subsequently, in a second step, at each spatial point, a state variable model using the reconstructed total concentration of the fluorophore in tissue as the measurement is used to estimate the states (compartment-concentrations) and the pharmacokinetic parameters. The estimation process uses an extended Kalman filter (EKF) framework.

Alternatively, in a one-step method, Alacam and Yazici⁹ and Wang et al.¹⁰ demonstrate that the use of log-intensity measurements to directly reconstruct spatially resolved compartment-concentrations and pharmacokinetic rates (rather than going via the pointwise total concentrations) offers potentially better reconstructions. Exhaustive reconstruction studies have been carried out by Alacam et al.⁹ comparing the performance of “one-step” linear as well as nonlinear inversions, with the “two-step”²¹ linearly modeled approach for diffusion-model studies in an EKF framework.

Dual-mode (x-ray CT with FOT) dynamic FOT pharmacokinetic reconstructions have also been approached in one and two-step least-squares based inversions, with linear Born-approximation propagation models using compartment-model derived biexponential temporal solutions, with structural priors obtained from x-ray CT.¹³^,²³^–²⁵ The structural priors better specify tissue optical properties corresponding to the organs housing various regions in the image, as well as regularizing the reconstructions.²⁶ In such a dual-mode setting,¹² a Karhunen–Loeve transform (KLT) based reconstruction is first carried out with a linear Born-type measurement model for the time-varying concentrations, and then a biexponential temporal model is used to obtain corresponding pharmacokinetic parameters.

Shape-based tomographic reconstructions are gaining importance²⁷^–³⁵ for reducing search space dimensions and thus enhancing computational tractability. A B-spline parametrization is used to represent absorption distribution in a two-dimensional (2-D) diffuse photon density wave model, which is solved using a greedy type optimization approach.³⁶ An ellipsoid representation of the absorption anomalies is used in an optical tomography (OT) problem setting using the Gauss–Newton (GN) method with line search.²⁸ Spherical harmonic parametrization is used to represent diffusion and absorption coefficients in a three-dimensional (3-D) OT problem, which is then solved using a line search-based GN scheme.²⁹ In an implicit parametrized level-set framework, Aghasi et al.³⁷ solve the diffusion optical tomography problem with radial basis function (RBF) based parametric level sets (PALS). A Hermite interpolation-based RBF representation is used in a 2-D FOT problem by Naik et al.³⁸ to represent the fluorophore absorption coefficient at the excitation wavelength; the FOT problem is then solved in pointwise and shape-based frameworks using Levenberg–Marquardt/GN methods.

A dynamic optical tomography problem in straight path ray-tomography was also solved³⁹ using RBF level-set object-boundary representations and GN filter based estimation of the dynamic shape and optical parameters. Nonparameterized pointwise-specified level-sets are used to implicitly represent the boundary of time varying fluorescence yield⁴⁰ in fluorescence molecular tomography. The pointwise level-set function and the piecewise constant values of the yield for the different time instants and projection angles are reconstructed using a gradient descent method.⁴⁰ It should be noted that a state variable model is not used in the abovementioned work to relate the fluorescent yield at different time instants.

Considering that we require spatially resolved images of various pharmacokinetic rates and volume fractions, with the objective of making the compartment model-based dynamic pharmacokinetic reconstruction problem more computationally tractable, we have proposed an RBF level-set parameterized shape-based tomographic inversion scheme using a regularized trust region based GN filter in a diffusion-approximation modeled FOT setting. Preliminary results from such a formulation in a Levenberg–Marquardt setting have been recently presented.⁴¹ In our pharmacokinetic tomographic settings, the (static) boundary of the tumor is reconstructed along with the dynamic concentrations within as well as the state-variable model’s pharmacokinetic parameters and volume fractions. Decay of the concentrations is assured by the structure of the state ordinary differential equation (ODE)-model’s coefficient-matrix. Hence, to ensure proper time-decay of concentrations, we directly solve for the pharmacokinetic parameters instead of going via their exponential propagator matrix components (designated as “interim kinetic parameters” in Ref. 10). Suitable error metrics have been then defined, and detailed numerical studies for near/sub-cm tumor-mimicking phantoms for two kinds of cancer and various data-SNRs have been presented and quantified with respect to the metrics. We see that our scheme yields a good localization of the test objects and reasonable estimates of the pharmacokinetic rates and concentration profiles.

The present paper differs from our recently presented work⁴¹ in that (a) the detailed formulation and derivation of the GN-filter Jacobians have been given in the present manuscript, (b) the details of the trust-region based regularized GN filter proposed have been given, and (c) error metrics have now been defined, and detailed numerical studies have now been included for tumor-mimicking phantoms of two different kinds of cancer [invasive ductal carcinoma (IDC) and adenocarcinoma (AC)] for various noise levels in the data and, importantly, the results have been well quantified by the error-metrics.

In Sec. 2, we discuss the compartment analysis of pharmacokinetics and shape representation of the states and parameters. Section 3 presents the proposed trust-region based iteratively regularized GN filter for the dynamic reconstruction as well as the Frechet derivatives of the measured log-intensity of the emission fluence with respect to pharmacokinetic and shape parameters. Section 4 contains the numerical studies on test cases of typical tumor mimicking phantoms, followed by the summary in Sec. 5. The Appendix gives the expressions for the interim kinetic parameters in terms of the pharmacokinetic parameters.

2. Problem Definition

2.1.

State Variable Model for Pharmacokinetic Compartment Analysis

Compartment modelling of pharmacokinetics is used in imaging studies⁹^,¹⁹^,²¹^,⁴²^–⁴⁴ for cancer detection. In compartment analysis, the region of interest is divided into virtual compartments or volumes, where the fluorophore concentration reaches rapid equilibrium upon injection.¹⁶^,⁴⁵ Indocyanine green (ICG) is an optical contrast agent, widely used for cancer detection studies.⁴⁶^–⁴⁸ A two compartment model (schematic in Fig. 1) has been reported to be suitable for describing ICG pharmacokinetics.⁴⁴ ICG administered intravenously into the blood stream binds to plasma proteins (albumin) and acts as a macromolecular agent.⁴³^,⁴⁶ Due to high leaky vasculature in the tumor region,⁴³^,⁴⁹ the macromolecule leaks into cancerous tissue, absorbs the incident light at the excitation wavelength, and emits (fluorescent) light at a longer wavelength, hence acting as a contrast agent for identifying tumors. The pharmacokinetic rates between the (blood) plasma compartment and the (tissue) extracellular and extravascular space compartment (EES) are higher in the tumor region due to the leaky nature of blood vessels there.

Fig. 1

Block diagram of two compartment model.

The fluorophore’s volume of distribution of a compartment is the apparent volume into which a given mass of fluorophore needs to be diluted to give the observed concentration.⁵⁰ As the contrast agent is leaked into the EES compartment, the apparent volume of distribution of ICG is more in the tumor region than in the healthy region. Thus, due to angiogenesis, defining the volume fractions of a compartment $(v_{e / p})$ as the ratio of its volume of distribution $(V_{e / p})$ and the total volume of distribution $(V = V_{e} + V_{p})$ , we see that the volume fraction of ICG in both compartments is greater in the tumor region.⁴³^,⁴⁴^,⁵¹ In the tumor regions, $v_{e}$ is in the range of 0.2 to 0.5,⁵² whereas $v_{p}$ is in the range of 0.013 to 0.067.⁴³^,⁴⁴

Let $C_{p} (μ M)$ and $C_{e} (μ M)$ be the concentration of ICG in the plasma compartment and EES compartment, respectively, $k_{p e} (s^{- 1})$ [respectively, $k_{p e} (s^{- 1})$ ] be the transfer rate of ICG from the plasma compartment to the EES compartment (respectively, the transfer rate of ICG from the EES compartment to the plasma compartment), and $k_{e l m} (s^{- 1})$ be the transfer rate at which ICG is eliminated from the region of interest.

The change in concentration of ICG in each compartment is described by the coupled ODE:⁹

Eq. (1)

\dot{C} (\vec{r}, t) = K (k_{p e} (\vec{r}), k_{e p} (\vec{r}), k_{e l m} (\vec{r})) C (\vec{r}, t),

where ‘

\vec{r}

’ denotes spatial coordinates and

C (\vec{r}, t) = [\begin{matrix} C_{e} (\vec{r}, t) \\ C_{p} (\vec{r}, t) \end{matrix}]

;

K (k_{p e} (\vec{r}), k_{e p} (\vec{r}), k_{e l m} (\vec{r})) = [\begin{matrix} - k_{e p} (\vec{r}) & k_{p e} (\vec{r}) \\ k_{e p} (\vec{r}) & - (k_{p e} (\vec{r}) + k_{e l m} (\vec{r})) \end{matrix}]

.

The corresponding discrete time-state model corresponding to time instants $t_{j}$ and $t_{j + 1}$ (indexed as $j$ and $j + 1$ , respectively) for Eq. (1) is given by⁹^,⁵³

Eq. (2)

C (\vec{r}, j + 1) = T (τ_{11} (\vec{r}), τ_{12} (\vec{r}), τ_{21} (\vec{r}), τ_{22} (\vec{r})) C (\vec{r}, j),

where

Eq. (3)

T \equiv e^{K t_{s}} \equiv [\begin{matrix} τ_{11} (\vec{r}) & τ_{12} (\vec{r}) \\ τ_{21} (\vec{r}) & τ_{22} (\vec{r}) \end{matrix}],

and

t_{s} = t_{j + 1} - t_{j}

is the sampling interval.

Denoting excitation and emission related quantities by subscripts “ $x$ ” and “ $m$ ”, respectively, the frequency domain governing equations, which describe light propagation in tissue, are given by⁵⁴

Eq. (4)

\begin{array}{l} - \nabla \cdot (D_{x} \nabla ϕ_{x}) + k_{x} ϕ_{x} = S_{x} \\ - \nabla \cdot (D_{m} \nabla ϕ_{m}) + k_{m} ϕ_{m} = β ϕ_{x} \end{array} in Ω,

being subject to the Robin boundary conditions:

Eq. (5)

\begin{array}{l} \vec{n} \cdot (D_{x} \nabla ϕ_{x}) + b_{x} ϕ_{x} = 0 \\ \vec{n} \cdot (D_{m} \nabla ϕ_{m}) + b_{m} ϕ_{m} = 0 \end{array} on \partial Ω,

where

Eq. (6)

D_{x / m} = \frac{1}{3 (μ_{a (x / m) i} + μ_{a (x / m) f} + μ_{s (x / m)}^{'})}, β = \frac{ϕ_{q} μ_{a x f}}{1 - i ω τ_{l}}, b_{x / m} = \frac{1 - R_{(x / m)}}{2 (1 + R_{(x / m)})}, k_{x / m} = \frac{i ω}{c} + μ_{a (x / m) i} + μ_{a (x / m) f},

where “

x / m

” stands for either “

x

” (excitation) or “

m

” (emission),

\vec{n}

is the outward normal to the boundary,

S_{x} (W / {cm}^{2})

is the excitation source with modulation frequency

ω (rad / s)

,

ϕ_{x} (W / {cm}^{2})

is the excitation fluence,

ϕ_{m} (W / {cm}^{2})

is the emission fluence,

D_{x / m}

is the diffusion coefficient at the excitation/emission wavelength,

k_{x / m}

is a decay coefficient,

μ_{a (x / m) i}

and

μ_{a (x / m) f}

are the absorption coefficients due to intrinsic chromophores and extrinsic fluorophores, respectively,

μ_{s (x / m)}^{'}

is the respective reduced scattering coefficient (all in

{cm}^{- 1}

) at the two wavelengths,

β

and

ϕ_{q}

are the unitless emission source coefficient and fluorescence quantum efficiency, respectively,

τ_{l}

is the fluorescence lifetime(in s),

c

is the speed of light in the medium (cm/s),

i = \sqrt{- 1}

,

b_{x / m}

are Robin boundary coefficients, and

R_{x / m}

are the reflection coefficients. We use a frequency domain modeling of the FOT process because of the inherent advantage of such systems in time-sampling based applications especially at reasonably good data SNRs.

The measurements are the complex log-intensity (defined below) at the detector locations at the excitation and emission wavelengths in general. In the present work, we focus on obtaining the absorption coefficient of the tissue at the excitation wavelength $(μ_{a x f})$ from the measured complex log-intensity at the emission wavelength.⁵⁴ We assume an a priori known linear relationship between $μ_{a x f}$ and $μ_{a m f}$ .⁵⁴ The measured log-intensity is as follows:

Eq. (7)

\log (ϕ_{m} [{{\vec{r}}_{d}}]) \equiv \log | ϕ_{m} ({\vec{r}}_{d}) | + i ν ({\vec{r}}_{d}),

where

ϕ_{m} \equiv | ϕ_{m} ({\vec{r}}_{d}) | e^{i ν}

and

{\vec{r}}_{d}

denotes a detector location. Hence, we can formally express the discrete-domain measurement equation at a time instant

t

as follows:

Eq. (8)

y ({\vec{r}}_{d}, t) \equiv \log ϕ_{m} [{{\vec{r}}_{d}}, t] \equiv G ({\underline{μ}}_{a x f} [\vec{r}, t]),

where

y ({\vec{r}}_{d}, t)

represents the vector of measurements at time

t

,

{\underline{μ}}_{a x f} (t)

is the vector of unknown absorption coefficients on the computational grid at time

t

, and

G (\cdot)

represents the measurement operator; in our case,

G (\cdot)

is evaluated using the finite-element method for the solution of the governing fluorescence diffusion model in Eqs. (4) and (5).⁵⁴

The relation between the total fluorophore concentration in the region of interest and the absorption coefficient is given by⁹

Eq. (9)

μ_{a (x / m) f} (\vec{r}, t) = \ln 10 \cdot ϵ_{(x / m)} \cdot C (\vec{r}, t),

where

ϵ_{x / m}

is the fluorophore extinction coefficient and

C (\vec{r}, t)

is the total fluorophore concentration in the tissue and is given by⁹

Eq. (10)

C (\vec{r}, t) = v_{p} (\vec{r}) C_{p} (\vec{r}, t) + v_{e} (\vec{r}) C_{e} (\vec{r}, t) .

2.2.

Level-Set Representation and State Variable Model

We consider a level-set based representation of a pharmacokinetic parameter $k_{ξ}$ , with $ξ \in {e p, p e, e l m}$ , as follows:

Eq. (11)

k_{ξ} (\vec{r}) = k_{ξ}^{i} (\vec{r}) H [ϕ_{γ} (\vec{r})] + k_{ξ}^{o} (\vec{r}) [1 - H (ϕ_{γ} (\vec{r}))],

where

k_{ξ}^{i} (\cdot)

and

k_{ξ}^{o} (\cdot)

represent

k_{ξ} (\cdot)

values inside and outside the region of interest, respectively. In our work, we consider

k_{ξ}^{i}

and

k_{ξ}^{o}

to be respective constants inside and outside the tumor region.

H

represents the Heaviside function, and

ϕ_{γ} (\vec{r})

is a function whose zero level-set represents the boundary of the tumor, and whose value is positive (negative, respectively) inside (outside, respectively) the tumor region.

γ

represents a set of parameters that define this level-set. In our work, we are using a RBF representation of the object boundary with compactly supported PALS.³⁷

The compactly supported RBF level-set function can be written as³⁷ follows:

Eq. (12)

ϕ_{γ} (\vec{r}) \equiv ϕ (\vec{r}, \underset{γ}{\underset{⏟}{[α, ζ, χ]}}) = \sum_{l = 1}^{m} α_{l} ψ ({‖ ζ_{l} (\vec{r} - χ_{l}) ‖}^{†}),

Eq. (13)

\tilde{r} \equiv {‖ \vec{r} ‖}^{†} = \sqrt{{‖ \vec{r} ‖}^{2} + v^{2}},

where

v

is a small real number,

ψ

is a compactly supported RBF,⁵⁵

m

is the number of RBFs used,

α_{l}

is the weighting factor,

ζ_{l}

is the dilation factor, and

χ_{l}

denotes the RBF center coordinates.

In the present work, we use as basis functions a $C^{2}$ polynomial of degree 5:³⁷^,⁵⁵

Eq. (14)

ψ (\tilde{r}) = {(1 - \tilde{r})}_{+}^{4} (4 \tilde{r} + 1),

where the

{(\cdot)}_{+}

denotes a cutoff function⁵⁵ defined as follows:

Eq. (15)

{(x)}_{+} = x \cdot H (x),

where

H (x)

is the Heaviside function. We also use a mollified Heaviside function in our work:⁵⁶

Eq. (16)

H_{ϵ} (ϕ) = {\begin{cases} 1 & if ϕ > ϵ_{w}, \\ 0 & if ϕ < - ϵ_{w}, \\ \frac{1}{2} + \frac{ϕ}{2 ϵ_{w}} + \frac{1}{2 π} \sin (\frac{π ϕ}{ϵ_{w}}) & if | ϕ | \leq ϵ_{w}, \end{cases}

where

ϵ_{w}

is the half-width of the transition region of the Heaviside. Concentrations in different compartments

C_{e}

,

C_{p}

, which are dependent on pharmacokinetic rates, are similarly assumed to be piecewise constant, and they can be expressed as follows:

Eq. (17)

C_{(e / p)} (\vec{r}, j) = C_{(e / p)}^{i} (j) H_{ϵ} [ϕ_{γ} (\vec{r})] + C_{(e / p)}^{o} (j) [1 - H_{ϵ} (ϕ_{γ} (\vec{r}))] .

The volume fractions $v_{e}$ and $v_{p}$ being different in healthy and tumor regions can also be similarly expressed via their inside/outside values denoted as $v_{e}^{i / o}$ and $v_{p}^{i / o}$ , respectively, as follows:

Eq. (18)

v_{(e / p)} (\vec{r}) = v_{(e / p)}^{i} (\vec{r}) H_{ϵ} [ϕ_{γ} (\vec{r})] + v_{(e / p)}^{o} (\vec{r}) [1 - H_{ϵ} (ϕ_{γ} (\vec{r}))] .

Substituting the level-set representation of concentrations and volume fractions in the above equation (for relation between $μ_{a x f}$ and $C$ ), we obtain the following:

Eq. (19)

μ (\vec{r}, j) = \ln 10 \cdot ϵ \cdot [(C_{e} (\vec{r}, j) v_{e} (\vec{r}) + C_{p} (\vec{r}, j) v_{p} (\vec{r}))] = \ln 10 \cdot ϵ \cdot [(C_{e} {(j)}^{i} v_{e}^{i} + C_{p}^{i} (j) v_{p}^{i}) H_{ϵ} (ϕ_{γ} (\vec{r})) + (C_{e}^{o} (j) v_{e}^{o} + C_{p}^{o} (j) v_{p}^{o}) (1 - H_{ϵ} (ϕ_{γ} (\vec{r})))],

where the subscripts on

μ

and

ϵ

have been omitted for ease of notation.

Using the level-set representation of pharmacokinetic rates and concentrations, the coupled ODE [Eq. (1)] is rewritten as follows:

Eq. (20)

[\begin{matrix} {\dot{C}}_{e}^{i} (t) \\ {\dot{C}}_{p}^{i} (t) \end{matrix}] H_{ϵ} (ϕ_{γ} (\vec{r})) + [\begin{matrix} {\dot{C}}_{e}^{o} (t) \\ {\dot{C}}_{p}^{o} (t) \end{matrix}] [1 - H_{ϵ} (ϕ_{γ} (\vec{r}))] = [\begin{matrix} - k_{e p}^{i} & k_{p e}^{i} \\ k_{e p}^{i} & - (k_{p e}^{i} + k_{e l m}^{i}) \end{matrix}] [\begin{matrix} C_{e}^{i} (t) \\ C_{p}^{i} (t) \end{matrix}] H_{ϵ} [ϕ_{γ} (\vec{r})] + [\begin{matrix} - k_{e p}^{o} & k_{p e}^{o} \\ k_{e p}^{o} & - (k_{p e}^{o} + k_{e l m}^{o}) \end{matrix}] [\begin{matrix} C_{e}^{o} (t) \\ C_{p}^{o} (t) \end{matrix}] [1 - H_{ϵ} (ϕ_{γ} (\vec{r}))] .

The time-discretized version for the above equation can thus be written as follows:

Eq. (21)

[\begin{matrix} C_{e}^{i} (j + 1) \\ C_{p}^{i} (j + 1) \\ C_{e}^{o} (j + 1) \\ C_{p}^{o} (j + 1) \end{matrix}] = [\begin{matrix} τ_{11}^{i} & τ_{12}^{i} & 0 & 0 \\ τ_{21}^{i} & τ_{22}^{i} & 0 & 0 \\ 0 & 0 & τ_{11}^{o} & τ_{12}^{o} \\ 0 & 0 & τ_{21}^{o} & τ_{22}^{o} \end{matrix}] [\begin{matrix} C_{e}^{i} (j) \\ C_{p}^{i} (j) \\ C_{e}^{o} (j) \\ C_{p}^{o} (j) \end{matrix}],

where the relations between the interim kinetic parameters, the

τ

’s, and the pharmacokinetic parameters, the

k

’s, arise from Eq. (3) and are given in the Appendix.

This state (concentration) and parameter (pharmacokinetic as well as shape parameters) estimation problem can be solved using either stochastic estimation schemes, such as EKF and its variants,⁹^,¹⁰^,²²^,⁵⁷ or deterministic schemes, such as the GN filter.³⁹^,⁵⁸ In our work, we propose an iteratively regularized deterministic GN-filter in a trust-region setting for our reconstructions. We define the state vector at time $j$ as follows:

Eq. (22)

Θ_{j} \equiv {\underset{C}{\underset{⏟}{C_{e}^{i} (j), C_{p}^{i} (j), C_{e}^{o} (j), C_{p}^{o} (j)}}, \underset{k}{\underset{⏟}{k_{p e}^{i}, k_{e p}^{i}, k_{e l m}^{i}, k_{p e}^{o}, k_{e p}^{o}, k_{e l m}^{o}}}, \underset{v}{\underset{⏟}{v_{e}^{i}, v_{e}^{o}, v_{p}^{i}, v_{p}^{o}}}, \underset{γ}{\underset{⏟}{α, β, χ_{1}, χ_{2}}}} .

Assuming our state equation is exact, we would need to explicitly reconstruct only $Θ_{0}$ . We can rewrite the state equation as follows:

Eq. (23)

[\begin{matrix} C_{e}^{i} (j + 1) \\ C_{p}^{i} (j + 1) \\ C_{e}^{o} (j + 1) \\ C_{p}^{o} (j + 1) \\ k \\ v \\ γ \end{matrix}] = [\begin{matrix} τ_{11}^{i} & τ_{12}^{i} & 0 & 0 & 0 & 0 & 0 \\ τ_{21}^{i} & τ_{22}^{i} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & τ_{11}^{o} & τ_{12}^{o} & 0 & 0 & 0 \\ 0 & 0 & τ_{21}^{o} & τ_{22}^{o} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & I_{6} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & I_{4} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & I_{4 m} \end{matrix}] [\begin{matrix} C_{e}^{i} (j) \\ C_{p}^{i} (j) \\ C_{e}^{o} (j) \\ C_{p}^{o} (j) \\ k \\ v \\ γ \end{matrix}],

where

I_{D}

denotes identity matrix of size

D \times D (D \in {6,4, 4 m})

. The above equation can be written in simplified form as follows:

Eq. (24)

Θ_{j + 1} = A (Θ_{j}) \cdot Θ_{j} = f (Θ_{j}),

where

f

denotes the nonlinear state transition function and matrix

A (Θ_{j})

is

Eq. (25)

A (Θ_{j}) \equiv [\begin{matrix} T^{i} & 0 & 0 \\ 0 & T^{o} & 0 \\ 0 & 0 & I_{4 m + 10} \end{matrix}]; T^{i / o} = [\begin{matrix} τ_{11}^{i / o} (k) & τ_{12}^{i / o} (k) \\ τ_{21}^{i / o} (k) & τ_{22}^{i / o} (k) \end{matrix}] .

The discrete-time measurement equation at time $j$ can be formally written from Eq. (8) as follows:

Eq. (26)

y_{j} \equiv g_{j} (Θ_{j}) = g_{j} [f_{j - 1} (\dots f_{0} (Θ_{0}))] .

3. Reconstruction Scheme

3.1.

Gauss–Newton Filter Scheme

The GN filter solves the state variable model with the state Eq. (23) and measurement Eq. (26) by solving the following regularized nonlinear least squares problem using the GN method:³⁹^,⁵⁸

Eq. (27)

{\hat{Θ}}_{0} = \underset{Θ_{0}}{\arg \min} F (Θ_{0}) ≔ \frac{1}{2} {‖ (g (Θ_{0}) - y) ‖}^{2} + τ R (Θ_{0}),

where

τ

is the regularization parameter,

R (\cdot)

is the regularization functional, and

y

and

g (Θ_{0})

are the concatenated set of observed and model-predicted measurements, respectively.

$R (Θ_{0})$ , the regularization functional, is chosen here as a minimum-norm penalty:

Eq. (28)

R (Θ_{0}) = {‖ Θ_{0} - Θ_{c} ‖}^{2},

where

Θ_{c}

represents an a priori known constant vector. The nonlinear least squares problem can be solved by an iterative regularization scheme using either line search or trust region methodologies.³⁹^,⁵⁹^,⁶⁰ A regularized GN update

p_{Θ}

solves at the current iterate

Θ

:

Eq. (29)

{\hat{p}}_{Θ} = \underset{p_{Θ}}{\arg \min} {‖ \begin{matrix} J (Θ) p_{Θ} + r \\ \sqrt{τ} (Θ - Θ_{c} + p_{Θ}) \end{matrix} ‖}^{2},

where the Jacobian

J

and the residual

r

are given by

Eq. (30)

J = [\begin{matrix} J_{0} \\ ⋮ \\ J_{M - 1} \end{matrix}], r = [\begin{matrix} r_{0} \\ ⋮ \\ r_{M - 1} \end{matrix}],

where

M

denotes the number of time instants and Jacobian at time instant

j

is given by

Eq. (31)

J_{j - 1} = G_{j - 1} [Θ_{j - 1}] F_{j - 1} [Θ_{j - 1}] F_{0} [Θ_{0}],

where

G [\cdot]

and

F [\cdot]

are the Frechet derivatives of measurement and state transition functions, respectively. The residual at time instant

j

is given by

Eq. (32)

r_{j} = g_{j} (Θ_{j}) - y_{j} = g_{j} [f_{j - 1} (\dots f_{0} (Θ_{0}))] - y_{j} .

3.2.

Reconstruction Algorithm

To solve the above minimization problem, we propose a trust region based iteratively regularized GN filter algorithm. We first observe that for a linear residual, the GN method converges in a single iteration. An iteratively (Tikhonov) regularized scheme needs to solve a succession of nonlinear regularized least squares subproblems, each with different values of the regularization parameter, with the solution obtained for a given parameter used as the starting estimate for the next lower parameter-value.

Our computational experience shows that we need not solve each subproblem exactly; hence, assuming that a linear assumption to the residual is approximately valid, we shift to lower parameter values if we are close to a full Newton step for a “good” update,³² with the “goodness” of the step being decided upon by the actual reduction in the residual with respect to its predicted decrease in our currently used trust region framework explained below.

Now, each step of a GN scheme needs to solve the problem in Eq. (29), which is a least-squares version of the linear system $J_{a} p_{Θ} = - r_{a}$ , where we use the augmented Jacobian $J_{a} \equiv [\begin{matrix} J \\ \sqrt{τ} I \end{matrix}]$ and the augmented residual $r_{a} \equiv [\begin{matrix} r \\ \sqrt{τ} (Θ_{0} - Θ_{c}) \end{matrix}]$ .

The step $p_{Θ}$ obtained from a GN step is found using either line-search or trust-region approaches so that the overall algorithm exhibits global convergence. In our work, we choose to work with the trust region class of schemes, wherein we assume a quadratic model of the cost-function to locally hold in some ball (decided by a trust region radius) around the current estimate. The trust region’s radius is varied with iteration as the ratio of the actual to predicted cost-function decreases.

Given a trust-region radius, the trust region step satisfies the following:

Eq. (33)

Δ = {‖ {(J_{a}^{T} J_{a} + λ I)}^{- 1} J_{a}^{T} r_{a} ‖}_{2},

where we calculate the parameter

λ

using a suitable root-finding technique.⁶⁰^,⁶¹

To ensure proper scaling, the variables are scaled as $Θ = S \cdot \tilde{Θ}$ using a diagonal scaling matrix $S$ , with elements given by⁵⁹

Eq. (34)

S_{i i} = \frac{1}{\sqrt{\sum_{j} J_{i j}^{2} + τ}} .

Thus, from Eq. (33), we see that the scaled-domain counterparts of $J_{a}$ and $Δ$ , namely, ${\tilde{J}}_{a} = J_{a} \cdot S$ and $\tilde{Δ}$ (the trust region radius in the scaled domain), respectively, satisfy

Eq. (35)

\tilde{Δ} = {‖ {({\tilde{J}}_{a}^{T} {\tilde{J}}_{a} + λ I)}^{- 1} {\tilde{J}}_{a}^{T} r_{a} ‖}_{2} .

The update in the scaled domain ${\tilde{p}}_{Θ}$ is then defined to satisfy the equation as follows:

Eq. (36)

[\begin{matrix} {\tilde{J}}_{a} \\ \sqrt{λ} I \end{matrix}] {\tilde{p}}_{Θ_{0}} = - [\begin{matrix} r_{a} \\ 0 \end{matrix}] .

The update in the original domain $p_{Θ}$ is thus given by

Eq. (37)

p_{Θ} = S {\tilde{p}}_{Θ} .

We then calculate the cost function $F$ and reduction ratio, $ρ = \frac{actual reduction}{predicted reduction}$ at the nominal estimate $Θ^{t}$ . The regularization parameter, $τ$ , is reduced if $ρ$ is above a threshold $ρ_{t h}$ with a minimum limit of $τ_{\min}$ . The trust region radius updating rule, values of $η_{1}$ , $η_{2}$ , and parameter $γ_{bad}$ are based on the practical algorithm in Conn et al.⁶² The algorithmic flow of reconstruction is shown in Algorithm 1.

Algorithm 1

Trust region-based iterative regularized GN filter.

1: Initialization:

Θ^{0}

,

Θ_{c} = Θ^{0}

,

Δ_{0}

,

η_{1} = 0.01

,

η_{2} = 0.9

,

τ_{0} = 0.8

,

i = 0

2: while $i < i_{max}$ do

3: Calculate

J^{i}

,

r^{i}

,

F (Θ^{i})

using

Θ^{i}

4: Calculate

λ

using Eq. (33)

5: Solve for

p_{Θ}

using Eqs. (36) and (37)

6:

Θ^{t} = Θ^{i} + p_{Θ}

7: Evaluate

F (Θ^{t})

and

ρ

8: if $ρ > η_{1}$ then

9: Accept update.

Θ^{i + 1} = Θ^{t}

;

10: if $ρ > ρ_{t h}$ then

11:

τ_{i + 1} = \max (τ_{i} / 3, τ_{\min})

;

12: end if

13: else

14:

Θ^{i + 1} = Θ^{i}

;

15: end if

16: if $ρ > η_{2}$ then

17:

Δ_{i + 1} = \max (2.5 \cdot ‖ {\tilde{p}}_{Θ} ‖, Δ_{i})

;

18: else if $ρ \geq η_{1} & ρ < η_{2}$ then

19:

Δ_{i + 1} = Δ_{i}

;

20: else if $ρ \geq 0 & ρ < η_{1}$ then

21:

Δ_{i + 1} = 0.25 \cdot ‖ {\tilde{p}}_{Θ} ‖

;

22: else if $ρ < 0$ then

23:

Δ_{i + 1} = \min (0.25 \cdot ‖ {\tilde{p}}_{Θ} ‖, \max (0.0625, γ_{b a d}) \cdot Δ_{i})

;

24: end if

25:

i = i + 1

;

26: end while

27: Choose the stopping iterate when

ϵ_{r e l} < t o l

or the data-residual stays stable or toggles across iterations.

The following relative measure³⁹ of the residual is used as a reflection of “how much” of the residual is left in the range of the augmented Jacobian and serves as a useful stopping criterion:

Eq. (38)

ϵ_{rel} = \frac{‖ P_{J_{a}} r_{a} ‖}{‖ r_{a} ‖},

where

P_{J_{a}}

is the orthogonal projection onto the range space of

J_{a}

.

3.3.

Frechet Derivative Calculation

3.3.1.

Frechet derivative of measurement function

The Frechet derivative of the measurement function [Eq. (8)] with respect to the parameter set at a given time instant $j$ is given by

Eq. (39)

G_{j} = [\begin{matrix} \underset{D \times 4}{\underset{⏟}{\frac{\partial y_{j}}{\partial C (j)}}} & \underset{D \times 6}{\begin{matrix} \underset{⏟}{\frac{\partial y_{j}}{\partial k}} \end{matrix}} & \underset{D \times 4}{\begin{matrix} \underset{⏟}{\frac{\partial y_{j}}{\partial v}} \end{matrix}} & \underset{D \times 4 m}{\begin{matrix} \underset{⏟}{\frac{\partial y_{j}}{\partial γ}} \end{matrix}} \end{matrix}] .

The sensitivity equation with respect to $s \in {C_{e}^{i} (j), C_{p}^{i} (j), C_{e}^{o} (j), C_{p}^{o} (j), k_{p e}^{i}, k_{e p}^{i}, k_{e l m}^{i}, k_{p e}^{o}, k_{e p}^{o}, k_{e l m}^{o}, v_{e}^{i}, v_{e}^{o}, v_{p}^{i}, v_{p}^{o}, α, β, χ_{1}, χ_{2}}$ is calculated using the chain rule as follows:

Eq. (40)

\frac{\partial y_{j} ({\vec{r}}_{d})}{\partial s} = \frac{1}{ϕ_{m} ({\vec{r}}_{d}, j)} [\frac{\partial ϕ_{m} ({\vec{r}}_{d}, j)}{\partial μ_{axf} (j)} \times \frac{\partial μ_{axf} (j)}{\partial s} + \frac{\partial ϕ_{m} ({\vec{r}}_{d}, j)}{\partial μ_{amf} (j)} \times \frac{\partial μ_{amf} (j)}{\partial s}] .

The derivative of the fluence $ϕ_{m}$ with respect to $μ_{axf}$ is given in Fedele et al.⁵⁴ To derive the sensitivity of $μ_{axf}$ and $μ_{amf}$ with respect to unknowns $Θ$ , consider Eq. (19) for the time-varying absorption coefficient:

Eq. (41)

μ (\vec{r}, j) = \ln 10 \cdot ϵ \cdot [(C_{e}^{i} (j) v_{e}^{i} + C_{p}^{i} (j) v_{p}^{i}) H_{ϵ} (ϕ (\vec{r})) + (C_{e}^{o} (j) v_{e}^{o} + C_{p}^{o} (j) v_{p}^{o}) (1 - H_{ϵ} (ϕ (\vec{r})))] .

The sensitivity of $μ (\vec{r}, j)$ with respect to ${C_{e}^{i} (j), C_{p}^{i} (j), C_{e}^{o} (j), C_{p}^{o} (j), v_{e}^{i}, v_{e}^{o}, v_{p}^{i}, v_{p}^{o}}$ is given below:

Eq. (42)

\frac{\partial μ (\vec{r}, j)}{\partial C_{e / p}^{i} (j)} = \ln 10 \cdot ϵ \cdot v_{e / p}^{i} H_{ϵ} [ϕ_{γ} (\vec{r})]; \frac{\partial μ (\vec{r}, j)}{\partial v_{e / p}^{i}} = \ln 10 \cdot ϵ \cdot C_{e / p}^{i} (j) H_{ϵ} [ϕ_{γ} (\vec{r})] .

Eq. (43)

\frac{\partial μ (\vec{r}, j)}{\partial C_{e / p}^{o} (j)} = \ln 10 \cdot ϵ \cdot v_{e / p}^{o} {1 - H_{ϵ} [ϕ_{γ} (\vec{r})]}; \frac{\partial μ (\vec{r}, j)}{\partial v_{e / p}^{o}} = \ln 10 \cdot ϵ \cdot C_{e / p}^{o} (j) {1 - H_{ϵ} [ϕ_{γ} (\vec{r})]} .

The variation of $μ (\vec{r}, j)$ with respect to $γ \in {α, ζ, χ}$ is given by

Eq. (44)

\frac{\partial μ (\vec{r}, j)}{\partial γ} = \frac{\partial μ (\vec{r}, j)}{\partial C_{e} (\vec{r}, j)} \frac{\partial C_{e} (\vec{r}, j)}{\partial γ} + \frac{\partial μ (\vec{r}, j)}{\partial C_{p} (\vec{r}, j)} \frac{\partial C_{p} (\vec{r}, j)}{\partial γ} + \frac{\partial μ (\vec{r}, j)}{\partial v_{e} (\vec{r})} \frac{\partial v_{e} (\vec{r})}{\partial γ} + \frac{\partial μ (\vec{r}, j)}{\partial v_{p} (\vec{r})} \frac{\partial v_{p} (\vec{r})}{\partial γ},

where [from Eq. (19)]

Eq. (45)

\frac{\partial μ (\vec{r}, j)}{\partial C_{e / p} (\vec{r}, j)} = \ln 10 \cdot ϵ \cdot v_{e / p} (\vec{r}), \frac{\partial μ (\vec{r}, j)}{\partial v_{e / p} (\vec{r})} = \ln 10 \cdot ϵ \cdot C_{e / p} (\vec{r}, j) .

Further, from the level-set representation of $C_{e / p}$ Eq. (17) and $v_{e / p}$ Eq. (18), we have

Eq. (46)

\frac{\partial C_{e / p} (\vec{r}, j)}{\partial γ} = [C_{e / p}^{i} (j) - C_{e / p}^{o} (j)] H_{ϵ}^{'} [ϕ_{γ} (\vec{r})] \frac{\partial ϕ}{\partial γ}, \frac{\partial v_{e / p} (\vec{r})}{\partial γ} = (v_{e / p}^{i} - v_{e / p}^{o}) H_{ϵ}^{'} [ϕ_{γ} (\vec{r})] \frac{\partial ϕ}{\partial γ},

where the sensitivities of the level-set,

ϕ

with repect to the shape parameters, are given by³⁷

Eq. (47)

\frac{\partial ϕ}{\partial α_{l}} = ψ ({‖ ζ_{j} (x - χ_{j}) ‖}^{†}),

Eq. (48)

\frac{\partial ϕ}{\partial ζ_{l}} = α_{j} ζ_{l} \frac{{‖ (x - χ_{l}) ‖}^{2}}{{‖ ζ_{l} (x - χ_{l}) ‖}^{†}} ψ^{'} ({‖ ζ_{j} (x - χ_{l}) ‖}^{†}),

Eq. (49)

\frac{\partial ϕ}{\partial χ_{l}^{k}} = - α_{l} ζ_{l}^{2} \frac{(x^{k} - χ_{l}^{k})}{{‖ ζ_{l} (x - χ_{l}) ‖}^{†}} ψ^{'} ({‖ ζ_{l} (x - χ_{l}) ‖}^{†}),

where

χ_{l}^{k}

is the

k

’th component of center

χ_{l}

.

3.3.2.

Frechet derivative of state transition function

The state transition equation is given by

Eq. (50)

[\begin{matrix} C_{e}^{i} (j + 1) \\ C_{p}^{i} (j + 1) \\ C_{e}^{o} (j + 1) \\ C_{p}^{o} (j + 1) \\ k (j + 1) \\ v (j + 1) \\ α (j + 1) \\ ζ (j + 1) \\ χ (j + 1) \end{matrix}] = [\begin{matrix} b_{1}^{i} (k, j) \\ b_{2}^{i} (k, j) \\ b_{1}^{o} (k, j) \\ b_{2}^{o} (k, j) \\ k (j) \\ v (j) \\ α (j) \\ ζ (j) \\ χ (j) \end{matrix}],

where

Eq. (51)

b_{1}^{i} (k, j) \equiv τ_{11}^{i} (k) C_{e}^{i} (j) + τ_{12}^{i} (k) C_{p}^{i} (j), b_{2}^{i} (k, j) \equiv τ_{21}^{i} (k) C_{e}^{i} (j) + τ_{22}^{i} (k) C_{p}^{i} (j), b_{1}^{o} (k, j) \equiv τ_{11}^{o} (k) C_{e}^{o} (j) + τ_{12}^{o} (k) C_{p}^{o} (j), b_{2}^{o} (k, j) \equiv τ_{21}^{o} (k) C_{e}^{o} (j) + τ_{22}^{o} (k) C_{p}^{o} (j) .

The variables $τ_{p q}^{i / o} (p, q = 1,2)$ are functions of parameters ${k_{e p}^{i / o}, k_{p e}^{i / o}, k_{e l m}^{i / o}}$ , whose expressions are given in the Appendix. The Frechet derivative, $F_{j}$ , of the state transition function at time instant $j$ is given by

F_{j} = [\begin{matrix} τ_{11}^{i} & τ_{21}^{i} & 0 & 0 & \frac{\partial b_{1}^{i}}{\partial k_{p e}^{i}} & \frac{\partial b_{1}^{i}}{\partial k_{e p}^{i}} & \frac{\partial b_{1}^{i}}{\partial k_{e l m}^{i}} & 0 & 0 & 0 & 0 \\ τ_{21}^{i} & τ_{22}^{i} & 0 & 0 & \frac{\partial b_{2}^{i}}{\partial k_{p e}^{i}} & \frac{\partial b_{2}^{i}}{\partial k_{e p}^{i}} & \frac{\partial b_{2}^{i}}{\partial k_{e l m}^{i}} & 0 & 0 & 0 & 0 \\ 0 & 0 & τ_{11}^{o} & τ_{21}^{o} & 0 & 0 & 0 & \frac{\partial b_{1}^{o}}{\partial k_{p e}^{o}} & \frac{\partial b_{1}^{o}}{\partial k_{e p}^{o}} & \frac{\partial b_{1}^{o}}{\partial k_{e l m}^{o}} & 0 \\ 0 & 0 & τ_{21}^{o} & τ_{22}^{o} & 0 & 0 & 0 & \frac{\partial b_{2}^{o}}{\partial k_{p e}^{o}} & \frac{\partial b_{2}^{o}}{\partial k_{e p}^{o}} & \frac{\partial b_{2}^{o}}{\partial k_{e l m}^{o}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & I_{4 m + 10} \end{matrix}] .

4. Numerical Studies

4.1.

Test-Case and Reconstruction Settings

A computational domain of size $4 \times 4 cm$ , with the origin as the center, is considered for our numerical test cases.⁶³ Eight detectors are placed symmetrically on each of the four sides of the domain, as shown in Fig. 2. Four collimated sources, each with strength 10 mW modulated at 100 MHz, are placed at the center of each side and modelled at the depth of one mean free path as in Ref. 64. The homogeneous optical properties of the tissue-mimicking phantom used in FOT with ICG as a fluorophore (at the excitation and emission wavelengths 785 and 830 nm, respectively) are given by⁶⁵ $μ_{axi} = 0.031 {cm}^{- 1}$ , $μ_{ami} = 0.00415 {cm}^{- 1}$ , $μ_{sx}^{'} = 10.95 {cm}^{- 1}$ , $μ_{sm}^{'} = 9.29 {cm}^{- 1}$ , $τ = 0.56 ns$ , $ϕ = 0.016$ , $R_{x, m} = 0.431$ , $ϵ_{x} = 130,000 M^{- 1} {cm}^{- 1}$ , $ϵ_{m} = 11,000 M^{- 1} {cm}^{- 1}$ .

Fig. 2

(a) Source and detector setting for the numerical study and (b) a typical fluence map obtained for the two object phantom for the IDC case (parameters in Table 2) with source on bottom face at $(0, - 2)$ .

The pharmacokinetic rates mentioned for IDC and AC²¹ are used to obtain synthetic measurement data. We assume $6.5 μ M$ concentration of fluorophore is injected⁴³ via bolus. At the first time instant in the data generation, we assume the fluorophore concentration in the plasma compartment to be 6.5 and $0 μ M$ in the EES compartment. Measurements are taken for 40 time instants with a sampling interval of 5 s. At each time instant, one source is on and measurements are taken from all the detectors (32 in our setting). Measurements used for reconstruction are the complex log intensity⁵⁹ of the fluence at the emission wavelength for all time instants (1280 in this setting).

The data are generated using a finer mesh discretized with 160,801 nodes containing 320,000 triangular elements. Reconstructions are performed using a coarser mesh discretized with 6561 nodes containing 12,800 triangular elements. Frechet derivatives evaluated using the method of adjoints are validated using the finite-difference method.

Numerical studies are done for each of the two cancer types (IDC and AC), for two phantoms: “T” denoting a two-object phantom (adjacently placed “smoothed corner square like” objects) with each having approximate extent of 0.5 cm in each direction with their boundaries separated by $\sim 1.5 cm$ and “B” being a single bean shaped phantom with lateral and longitudinal extents being $\sim 0.7$ and 1.3 cm, respectively. Data are generated for the two phantoms at three SNR levels, as given in Table 1. All the computations are performed in the Matlab^® 2016a programming environment.

Table 1

SNR of the synthetic data.

T phantom (IDC)	SNR (dB)	B phantom (IDC)	SNR (dB)
I-T1	38.99	I-B1	39.35
I-T2	33.2	I-B2	33.47
I-T3	29.77	I-B3	29.75
T Phantom (AC)	SNR (dB)	B phantom (AC)	SNR (dB)
A-T1	39.02	A-B1	39.29
A-T2	32.87	A-B2	33.4
A-T3	29.35	A-B3	29.83

Initial estimates for pharmacokinetic rates are taken in between healthy and tumor values. Initial fluorophore concentration in plasma compartment ( $C_{p}^{i}$ and $C_{p}^{o}$ ) is assumed to be $6.5 μ M$ and concentration of the fluorophore in EES compartment ( $C_{e}^{i}$ and $C_{e}^{o}$ ) is assumed to be $0 μ M$ . The algorithm is terminated when $ϵ_{rel} < tol$ or the data-residual remains stable or toggles across iterations.

The shape reconstructions of two-object and bean phantoms for IDC as well as AC under various noise conditions are shown in Figs. 3 and 4, respectively. The concentration curves (with respect to time) for both phantoms corresponding to regions inside and outside the tumor in IDC settings are shown in Fig. 5 for two-object phantom and Fig. 6 for bean phantom and in AC settings are shown in Fig. 7 for two-object phantom and Fig. 8 for bean phantom. In the figures, a time-index is defined as representing the sampling interval used (5 s in our case).

Fig. 3

Reconstruction of two object phantoms for IDC (a) dataset I-T1, (c) dataset I-T2, (e) dataset I-T3 and AC (b) dataset A-T1, (d) dataset A-T2, (f) dataset A-T3 tumor cases. Blue dotted line denotes the initial level-set, red dashed line denotes the shape of true object, and black solid line denotes the reconstructed shape.

Fig. 4

Reconstruction of bean shape phantom for IDC (a) dataset I-B1, (c) dataset I-B2, (e) dataset I-B3 and AC (b) dataset A-B1, (d) dataset A-B2, (f) dataset A-B3 tumor cases. Blue dotted line denotes the initial level-set, red dashed line denotes the shape of true object, and black solid line denotes the reconstructed shape.

Fig. 5

Concentration versus time plot for two object phantom for IDC tumor; 1 time-index = 5 s. (a) Schematic of phantom placed on top. Red denotes the decay of concentration in true phantom and blue denotes the decay in reconstructed phantom. (b, d, f) Left column shows $C e_{in}$ and $C p_{in}$ plots in the tumor region. (c, e, g) Right column shows $C e_{out}$ and $C p_{out}$ plots outside the tumor region. (b, c) First row corresponds to dataset I-T1, (d, e) second row to dataset I-T2, and (f, g) third row to dataset I-T3.

Fig. 6

Concentration versus time plot for bean shape object phantom for IDC tumor; 1 time-index = 5 s. (a) Schematic of phantom placed on top. Red denotes the decay of concentration in true phantom and blue denotes the decay in reconstructed phantom. (b, d, f) Left column shows $C e_{in}$ and $C p_{in}$ plots in the tumor region. (c, e, g) Right column shows $C e_{out}$ and $C p_{out}$ plots outside the tumor region. (a, b) First row corresponds to dataset I-B1, (c, d) second row to dataset I-B2, and (e, f) third row to dataset I-B3.

Fig. 7

Concentration versus time plot for two object phantom for AC tumor; 1 time-index = 5 s. (a) Schematic of phantom placed on top. Red denotes the decay of concentration in true phantom and blue denotes the decay in reconstructed phantom. (b, d, f) Left column shows $C e_{in}$ and $C p_{in}$ plots in the tumor region. (c, e, g) Right column shows $C e_{out}$ and $C p_{out}$ plots outside the tumor region. (b, c) First row corresponds to dataset A-T1, (d, e) second row to dataset A-T2, and (f, g) third row to dataset A-T3.

Fig. 8

Concentration versus time plot for bean shape object phantom for AC tumor; 1 time-index = 5 s. (a) Schematic of phantom placed on top. Red denotes the decay of concentration in true phantom and blue denotes the decay in reconstructed phantom. (b, d, f) Left column shows $C e_{in}$ and $C p_{in}$ plots in the tumor region. (c, e, g) Right column shows $C e_{out}$ and $C p_{out}$ plots outside the tumor region. (b, c) First row corresponds to dataset A-B1, (d, e) second row to dataset A-B2, and (f, g) third row to dataset A-B3.

Tables 2 and 3 show the reconstructed values of the pharmacokinetic parameters for IDC and AC tumor cases, respectively. To gauge the performance of the algorithm, across all test cases considered (i.e., along each row), we evaluate across all cases for each parameter (a) the normalized mean square error (row NMSE) of the reconstructions and (b) the maximal normalized error (MNE) in the reconstruction. The evaluations in (a) and (b) above are found in the last two columns, respectively, of Tables 2 and 3.

Table 2

Pharmacokinetic parameters for IDC test cases with average and maximal errors.

Reconstructed values
Parameter	True	I-T1	I-T2	I-T3	I-B1	I-B2	I-B3	Row NMSE	Maximal normalized error
$C_{e}^{i}$	0	0.01	0.01	0.01	$6 \times 10^{- 4}$	0.01	0.01	$8 \times 10^{- 5}$ a	0.01a
$C_{e}^{o}$	0	0	0.0029	0	0	0	0	$1 \times 10^{- 6}$ a	0.0029a
$C_{p}^{i}$	6.5	6.49	6.5	6.5	6.49	6.5	6.45	$8 \times 10^{- 6}$	0.0073
$C_{p}^{o}$	6.5	6.5	6.5	6.5	6.5	6.46	6.5	$5 \times 10^{- 6}$	0.0059
$k_{p e}^{i}$	0.0687	0.0896	0.0900	0.0900	0.0453	0.0803	0.0659	0.07	0.34
$k_{p e}^{o}$	0.0306	0.0282	0.0282	0.0265	0.0315	0.0225	0.0220	0.03	0.28
$k_{e p}^{i}$	0.0496	0.0436	0.0583	0.0470	0.0446	0.0586	0.0484	0.01	0.18
$k_{e p}^{o}$	0.0166	0.0142	0.0140	0.0206	0.0185	0.0153	0.0157	0.02	0.24
$k_{e l m}^{i}$	0.00449	0.0070	0.0047	0.0070	0.0041	0.0028	0.0042	0.13	0.55
$k_{elm}^{o}$	0.00446	0.0029	0.0047	0.0038	0.0041	0.0028	0.0040	0.05	0.38
$v_{e}^{i}$	0.3	0.4423	0.4390	0.4770	0.3556	0.3144	0.3385	0.14	0.58
$v_{e}^{o}$	0	0	$3.1 \times 10^{- 6}$	0	0	0	0	$1.5 \times 10^{- 12}$ a	$3.1 \times 10^{- 6}$ a
$v_{p}^{i}$	0.0600	0.0700	0.0498	0.0700	0.0669	0.0700	0.0700	0.02	0.16
$v_{p}^{o}$	0.0200	0.0190	0.0200	0.0160	0.0198	0.0167	0.0175	0.01	0.19

^aThe error is not normalized owing to true values being zero.

Table 3

Pharmacokinetic parameters for AC test cases with average and maximal errors.

Reconstructed values
Parameter	True	A-T1	A-T2	A-T3	A-B1	A-B2	A-B3	Row NMSE	Maximal normalized error
$C_{e}^{i}$	0	0.01	0.01	0.01	0.01	0.0017	0	$8.37 \times 10^{- 5}$ a	0.01a
$C_{e}^{o}$	0	0	0	0	0	0	0	0a	0a
$C_{p}^{i}$	6.5	6.5	6.5	6.5	6.5	6.49	6.5	$2 \times 10^{- 7}$	$1.3 \times 10^{- 3}$
$C_{p}^{o}$	6.5	6.5	6.5	6.5	6.49	6.5	6.5	$1 \times 10^{- 10}$	$3 \times 10^{- 5}$
$k_{p e}^{i}$	0.0292	0.0225	0.0199	0.0229	0.0223	0.0286	0.0216	0.05	0.31
$k_{p e}^{o}$	0.0114	0.0103	0.0136	0.0133	0.0102	0.0071	0.0075	0.05	0.37
$k_{e p}^{i}$	0.0158	0.0118	0.0136	0.0117	0.0100	0.0121	0.0100	0.08	0.36
$k_{e p}^{o}$	0.0065	0.0060	0.0066	0.0082	0.0069	0.0067	0.0057	0.02	0.26
$k_{elm}^{i}$	0.0043	0.0038	0.0025	0.0031	0.0039	0.0032	0.0033	0.06	0.41
$k_{elm}^{o}$	0.0035	0.0035	0.0025	0.0031	0.0039	0.0032	0.0033	0.02	0.28
$v_{e}^{i}$	0.2000	0.0794	0.0645	0.0891	0.1317	0.0901	0.0967	0.3	0.67
$v_{e}^{o}$	0	0	0	0	0	0	0	0a	0a
$v_{p}^{i}$	0.0400	0.0472	0.0291	0.0342	0.0428	0.0416	0.0377	0.02	0.27
$v_{p}^{o}$	0.0200	0.0191	0.0199	0.0187	0.0196	0.0170	0.0185	0.01	0.15

^aThe error is not normalized owing to true values being zero.

The contrast in the pharmacokinetic-rates and plasma-volume fractions indicates the presence of high-permeability and angiogenesis, respectively, in the tumor region. From our reconstructions, we see that we are able to obtain a good delineation between tumor and healthy tissue as well as the contrast in the pharmacokinetic rates and plasma-volume fraction. We observe that the reconstructions are reasonably stable for datasets of SNR about 30 dB, below which the quality and stability of reconstructions tend to go down.

4.2.

Quantification of Errors

To quantify the quality of our shape-based reconstructions, we use four error measures, namely, normalized error of area-parameter product across time instants, the distance of the centroid of the reconstructed object from the actual object, the Dice coefficient for the shape reconstructions, and NMSE for the pharmacokinetic rates.

In addition, we also use an NMSE for pointwise evaluated pharmacokinetic rate and volume fraction images to get an image-quality metric for our shape-based reconstructions. We map our shape and pharmacokinetic parameter reconstructions into pointwise values to compute these NMSEs. The spatial values are obtained using Eq. (11).

The area-parameter product error measure,³⁸ $E_{AP}$ , is defined across all the time instants at which measurements are taken as follows:

Eq. (52)

E_{A P} = ([\frac{\sum_{j = 1}^{M} | μ_{r e c}^{i} (j) A_{r e c} - μ_{a c}^{i} (j) A_{a c} |}{\sum_{j = 1}^{M} | μ_{a c}^{i} (j) A_{a c} |}] / M) \times 100 %,

where

A_{a c}

and

A_{r e c}

represent the area of actual and reconstructed object, respectively,

μ_{r e c}^{i} (j)

[respectively,

μ_{a c}^{i} (j)

] represents reconstructed (respectively, actual) fluorophore absorption coefficient Eq. (41) inside the tumor region at time instant

j

. This definition allows emphasis with respect to time instants with more significant product values.

The area of an object is given by

Eq. (53)

A_{object} = a_{element} \sum_{i, j} χ_{object} (x_{i}, y_{j}),

where

a_{element}

is the area of an element (a constant in our studies),

χ_{object} (\cdot)

is the characteristic function with respect to object support, and

(i, j)

represents the indices of centroid coordinates

x

and

y

of the discretized domain. The centroid coordinates of an object are given by

Eq. (54)

{\bar{x}}_{object} = \frac{\sum_{i, j} x_{i} χ_{object} (x_{i}, y_{j})}{A_{object}}, {\bar{y}}_{object} = \frac{\sum_{i, j} y_{j} χ_{object} (x_{i}, y_{j})}{A_{object}} .

The Euclidean distance between the centroids of reconstructed object and actual object is given by

Eq. (55)

E_{c} = \sqrt{{({\bar{x}}_{r e c} - {\bar{x}}_{a c})}^{2} + {({\bar{y}}_{r e c} - {\bar{y}}_{a c})}^{2}} .

The Dice coefficient⁶⁶ quantifies the localization and similarity of the shape reconstruction with the original shape. If $S$ denotes the set of nodes inside the reconstructed object and $H$ denotes the set of nodes inside the true object, the Dice coefficient is given by

Eq. (56)

D (S, H) = \frac{2 | S \cap H |}{| S | + | H |},

where

| S \cap H |

denotes the number of nodes present in

S

and also belongs to

H

. The Dice coefficient varies from 0 (indicating complete mismatch) to 1 (indicating accurate shape reconstruction).

The NMSE defined for a reconstructed quantity $X^{r}$ with respect to its actual value $X^{a}$ is defined as follows:

Eq. (57)

E_{NMSE} = \frac{{‖ X^{r} - X^{a} ‖}^{2}}{{‖ X^{a} ‖}^{2}} .

The error metrics for the reconstruction of two phantoms at different noise levels are tabulated in Table 4. We note that the NMSE for the pharmacokinetic rates evaluated in this table is based on the reconstructions of the concatenated vector $k$ , as shown in Eq. (22).

Table 4

Error measures for the shape reconstructions [for the two object case the centroid error is an ordered pair (a,b) corresponding to each object].

Phantom	EAP (%)	EC (cm)	D(S,H)	ENMSE(k)
I-T1	1.29	(0.04, 0.04)	0.84	0.058
I-T2	2	(0.14, 0.07)	0.80	0.064
I-T3	1.24	(0.04, 0.08)	0.81	0.059
I-B1	0.2	0.01	0.88	0.068
I-B2	$8.9 \times 10^{- 2}$	$2 \times 10^{- 3}$	0.89	0.034
I-B3	0.34	0.07	0.79	0.01
A-T1	0.47	(0.02,0.03)	0.59	0.047
A-T2	0.66	(0.12,0.05)	0.46	0.076
A-T3	0.97	(0.21,0.03)	0.5	0.049
A-B1	0.51	0.02	0.91	0.063
A-B2	0.59	0.012	0.83	0.026
A-B3	0.85	0.07	0.84	0.082

The error metrics evaluated further emphasize the good localization in general given by our approach for the small (near/sub-cm) phantoms in our study, in addition to a reasonable error in the reconstructed pharmacokinetic parameters and the area-parameter-product aspect.

In Table 5, to relate our shape-based results to pointwise error estimates with the purpose of checking image-quality acceptability with respect to the existing literature (to the best of our knowledge, the only paper that solves for the present pharmacokinetic parameters along with volume fractions in a “one-step” reconstruction considering a fully nonlinear FOT model is the work of Alacam et al.;⁹ they solve the pointwise problem with an EKF), we evaluate NMSE values [in the form of $20 \log (NMSE) dB$ ] for the mapped-pointwise reconstructed images of the pharmacokinetic parameters and volume fractions. Our obtained pointwise NMSEs for $k_{p e}$ ( $k_{in}$ in Alacam et al.⁹) range from $- 31$ to $- 16.9 dB$ and those for $k_{e p}$ ( $k_{out}$ in Alacam et al.⁹) range from $- 31.3$ to $- 20.94 dB$ across data SNR levels. The work of Alacam et al.⁹ reports NMSEs of $- 19.77$ and $- 18.49 dB$ for $k_{in}$ and $k_{out}$ , respectively, for noiseless data with their synthetic phantoms; they do not report any error values for their volume fractions. This shows that our reconstructions are well within accepted ranges for reconstruction quality.

Table 5

Values of 20 log(NMSE)dB for the spatial pharmacokinetic parameter and volume fraction values.

Phantom	E(kpe) dB	E(kep) dB	E(ve) dB	E(vp) dB	E(kelm) dB
I-T1	$- 31$	$- 24.22$	$- 17.03$	$- 63.82$	$- 18.55$
I-T2	$- 30.9$	$- 25.07$	$- 17.36$	$- 67.11$	$- 51.5$
I-T3	$- 26.9$	$- 22.01$	$- 10.25$	$- 49.57$	$- 31.11$
I-B1	$- 30$	$- 28.9$	$- 26.2$	$- 65$	$- 44.4$
I-B2	$- 22.6$	$- 28.57$	$- 38.53$	$- 59$	$- 16.8$
I-B3	$- 20.76$	$- 25.74$	$- 20.13$	$- 51$	$- 40.05$
A-T1	$- 28.26$	$- 30.11$	$- 13.7$	$- 60$	$- 61.49$
A-T2	$- 22.61$	$- 24.65$	$- 9.38$	$- 79.29$	$- 21.59$
A-T3	$- 23.5$	$- 20.94$	$- 10.38$	$- 73.3$	$- 35.18$
A-B1	$- 30.49$	$- 28.3$	$- 32.01$	$- 86.1$	$- 39.62$
A-B2	$- 16.92$	$- 31.3$	$- 19.31$	$- 62.73$	$- 39.74$
A-B3	$- 18.55$	$- 26.11$	$- 20.2$	$- 79.38$	$- 42.68$

5. Summary

In this work, we propose a shape-based dynamic tomographic reconstruction scheme for fluorescence-based pharmacokinetics using a regularized GN filter approach. The contribution of the present work is to represent spatially varying pharmacokinetic parameters, fluorophore concentrations, and volume fractions using compactly supported RBF-based level-set representations and derive the corresponding shape based Frechet derivatives for time-varying log intensity measurements of the optical signal. An iteratively regularized trust region based GN filter has been proposed to solve this reconstruction problem. It should be noted that we directly reconstruct pharmacokinetic rates, rather than the state equation propagator components (the interim kinetic parameters) as in some previous works.⁹^,¹⁰

Numerical studies with noisy synthetic data obtained from tumor mimicking numerical phantoms having near/sub-cm dimensions are presented, which validate the proposed scheme. The reconstructions demonstrate good localization and reasonable shape and optical parameter reconstructions, thus demonstrating the good potential of this methodology as an early cancer diagnostic. To obtain a pointwise reconstruction error measure, we mapped our shape and pharmacokinetic parameter reconstructions into pointwise values to compute pointwise-image-NMSEs; comparison of these pointwise-image-NMSEs with the errors reported in the literature for the pharmacokinetic rates shows that they are well within acceptable ranges.

The aim of our detailed computational study is to obtain a clear understanding of the numerical characteristics of our proposed algorithm before going to experimental data. Aspects related to application to in vivo settings, in addition to the three-dimensional modeling requirement, would be:

1. Characterization of the data-acquisition in terms of limited-data aspects, source-detector configurations (especially depending on the region to be interrogated as well as the object representation chosen) and detection numerical apertures, detector sensitivity, and temporal resolution possible to obtain sufficient data SNRs (for accurate reconstructions, since we observe that, in our work, the data SNRs would be needed to be above about 30 dB) would be needed while applying the algorithm in actual physical settings.
2. The present results are for scattering dominant media, where the diffusion approximation holds as the governing model of light propagation. However, for tissues that are absorption dominant over the wavelengths of use, we will have to go in for forward models, such as the full RTE⁶⁷ or approximations, such as the simplified spherical harmonics ( $S P_{n}$ ) ones.³⁸^,⁶⁸
3. The development of computationally efficient algorithms in 3-D and detailed reconstruction studies with respect to image representation and data-acquisition configurations, which is the subject of ongoing work.

6. Appendix

Expressions of $τ_{i j} i, j = {1,2}$ in terms of $k_{ξ}$ , with $ξ = {e p, p e, e l m}$ , are given by Eqs. (58)–(60) (for ease of notations, we have omitted the superscript $i / o$ on $τ$ and $k$ -related quantities):

Eq. (58)

τ_{12} = \frac{4 k_{e p} k_{p e} (e^{Γ_{2}} - e^{Γ_{1}})}{4 k_{e p} \sqrt{Ξ_{1}}}, τ_{21} = - \frac{k_{e p} e^{Γ_{1}} - k_{e p} e^{Γ_{2}}}{\sqrt{Ξ_{1}}},

Eq. (59)

τ_{11} = \frac{e^{Γ_{1}} \sqrt{Ξ_{1}} + e^{Γ_{2}} \sqrt{Ξ_{1}} - k_{e l m} e^{Γ_{1}} + k_{e l m} e^{Γ_{2}} + k_{e p} e^{Γ_{1}^{u}} - k_{e p} e^{Γ_{2}} - k_{p e} e^{Γ_{1}} + k_{p e} e^{Γ_{2}}}{2 \sqrt{Ξ_{1}}},

Eq. (60)

τ_{22} = \frac{e^{Γ_{1}} \sqrt{Ξ_{1}} + e^{Γ_{2}} \sqrt{Ξ_{1}} + k_{e l m} e^{Γ_{1}} - k_{e l m} e^{Γ_{2}} - k_{e p} e^{Γ_{1}^{u}} + k_{e p} e^{Γ_{2}} + k_{p e} e^{Γ_{1}} - k_{p e} e^{Γ_{2}}}{2 \sqrt{Ξ_{1}}},

where

Ξ_{1}

,

Ξ_{2}

,

Γ_{1}

, and

Γ_{2}

are given by the following expressions:

Eq. (61)

Γ_{1} = Ξ_{2} - \frac{t_{s} \sqrt{Ξ_{1}}}{2}; Γ_{2} = \frac{t_{s} \sqrt{Ξ_{1}}}{2} + Ξ_{2},

Eq. (62)

Ξ_{1} = k_{e l m}^{2} - 2 k_{e l m} k_{e p} + 2 k_{e l m} k_{p e} + k_{e p}^{2} + 2 k_{e p} k_{p e} + k_{p e}^{2},

Eq. (63)

Ξ_{2} = - \frac{k_{e l m} t_{s}}{2} - \frac{k_{e p} t_{s}}{2} - \frac{k_{p e} t_{s}}{2} .

The derivatives of $τ_{i j}$ with respect to $k_{ξ}$ have expressions that are too large to be included here.

Disclosures

The authors have no relevant financial interests in this article and no potential conflicts of interests to disclose.

References

1.

S. R. Meikle et al., “Parametric image reconstruction using spectral analysis of PET projection data,” Phys. Med. Biol., 43 (3), 651 –666 (1998). https://doi.org/10.1088/0031-9155/43/3/016 PHMBA7 0031-9155 Google Scholar

2.

B. Reutter, G. Gullberg and R. Huesman, “Kinetic parameter estimation from attenuated SPECT projection measurements,” IEEE Trans. Nucl. Sci., 45 (6), 3007 –3013 (1998). https://doi.org/10.1109/23.737657 IETNAE 0018-9499 Google Scholar

3.

R. H. Huesman et al., “Kinetic parameter estimation from SPECT cone-beam projection measurements,” Phys. Med. Biol., 43 (4), 973 –982 (1998). https://doi.org/10.1088/0031-9155/43/4/024 PHMBA7 0031-9155 Google Scholar

4.

O. Langer et al., “Combined PET and microdialysis for in vivo assessment of intracellular drug pharmacokinetics in humans,” J. Nucl. Med., 46 (11), 1835 –1841 (2005). JNMEAQ 0161-5505 Google Scholar

5.

M. E. Kamasak et al., “Direct reconstruction of kinetic parameter images from dynamic PET data,” IEEE Trans. Med. Imaging, 24 (5), 636 –650 (2005). https://doi.org/10.1109/TMI.2005.845317 ITMID4 0278-0062 Google Scholar

6.

G. Wang and J. Qi, “Generalized algorithms for direct reconstruction of parametric images from dynamic pet data,” IEEE Trans. Med. Imaging, 28 (11), 1717 –1726 (2009). https://doi.org/10.1109/TMI.2009.2021851 ITMID4 0278-0062 Google Scholar

7.

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

8.

A. Godavarty et al., “Fluorescence-enhanced optical imaging in large tissue volumes using a gain-modulated ICCD camera,” Phys. Med. Biol., 48 (12), 1701 –1720 (2003). https://doi.org/10.1088/0031-9155/48/12/303 PHMBA7 0031-9155 Google Scholar

9.

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

10.

X. Wang et al., “Direct reconstruction in CT-analogous pharmacokinetic diffuse fluorescence tomography: two-dimensional simulative and experimental validations,” J. Biomed. Opt., 21 (4), 046007 (2016). https://doi.org/10.1117/1.JBO.21.4.046007 JBOPFO 1083-3668 Google Scholar

11.

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

12.

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

13.

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

14.

X. Zhang et al., “Fast reconstruction of fluorophore concentration variation based on the derivation of the diffusion equation,” J. Opt. Soc. Am. A, 32 (11), 1993 –2001 (2015). https://doi.org/10.1364/JOSAA.32.001993 JOAOD6 0740-3232 Google Scholar

15.

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

16.

Y. Kwon, Handbook of Essential Pharmacokinetics, Pharmacodynamics and Drug Metabolism for Industrial Scientists, Springer Science & Business Media, New York (2001). Google Scholar

17.

S. E. Rosenbaum, Basic Pharmacokinetics and Pharmacodynamics: An Integrated Textbook and Computer Simulations, John Wiley & Sons, Hoboken, New Jersey (2016). Google Scholar

18.

H. Sun et al., “Imaging the pharmacokinetics of [F-18] FAU in patients with tumors: PET studies,” Cancer Chemother. Pharmacol., 57 (3), 343 –348 (2006). https://doi.org/10.1007/s00280-005-0037-0 CCPHDZ 0344-5704 Google Scholar

19.

M. Gurfinkel et al., “Pharmacokinetics of ICG and HPPH-car for the detection of normal and tumor tissue using fluorescence, near-infrared reflectance imaging: a case study,” Photochem. Photobiol., 72 (1), 94 –102 (2000). https://doi.org/10.1562/0031-8655(2000)072<0094:POIAHC>2.0.CO;2 PHCBAP 0031-8655 Google Scholar

20.

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 (7), 1357 –1368 (2005). https://doi.org/10.1364/JOSAA.22.001357 JOAOD6 0740-3232 Google Scholar

21.

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

22.

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

23.

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

24.

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

25.

G. Zhang et al., “Imaging of pharmacokinetic rates of indocyanine green in mouse liver with a hybrid fluorescence molecular tomography/x-ray computed tomography system,” J. Biomed. Opt., 18 (4), 040505 (2013). https://doi.org/10.1117/1.JBO.18.4.040505 JBOPFO 1083-3668 Google Scholar

26.

S. C. Davis et al., “Image-guided diffuse optical fluorescence tomography implemented with Laplacian-type regularization,” Opt. Express, 15 (7), 4066 –4082 (2007). https://doi.org/10.1364/OE.15.004066 OPEXFF 1094-4087 Google Scholar

27.

F. Santosa, “A level-set approach for inverse problems involving obstacles,” ESAIM: Control Optim. Calculus Var., 1 17 –33 (1996). https://doi.org/10.1051/cocv:1996101 Google Scholar

28.

M. E. Kilmer et al., “Three-dimensional shape-based imaging of absorption perturbation for diffuse optical tomography,” Appl. Opt., 42 (16), 3129 –3144 (2003). https://doi.org/10.1364/AO.42.003129 APOPAI 0003-6935 Google Scholar

29.

A. D. Zacharopoulos et al., “Three-dimensional reconstruction of shape and piecewise constant region values for optical tomography using spherical harmonic parametrization and a boundary element method,” Inverse Probl., 22 (5), 1509 –1532 (2006). https://doi.org/10.1088/0266-5611/22/5/001 INPEEY 0266-5611 Google Scholar

30.

O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Probl., 22 (4), R67 (2006). https://doi.org/10.1088/0266-5611/22/4/R01 INPEEY 0266-5611 Google Scholar

31.

S. R. Arridge et al., “Parameter and structure reconstruction in optical tomography,” J. Phys. Conf. Ser., 135 (1), 012001 (2008). https://doi.org/10.1088/1742-6596/135/1/012001 JPCSDZ 1742-6588 Google Scholar

32.

N. Naik et al., “A nonlinear iterative reconstruction and analysis approach to shape-based approximate electromagnetic tomography,” IEEE Trans. Geosci. Remote Sens., 46 (5), 1558 –1574 (2008). https://doi.org/10.1109/TGRS.2008.916077 IGRSD2 0196-2892 Google Scholar

33.

O. Dorn and D. Lesselier, “Level set methods for inverse scattering - some recent developments,” Inverse Probl., 25 (12), 125001 (2009). https://doi.org/10.1088/0266-5611/25/12/125001 INPEEY 0266-5611 Google Scholar

34.

F. Larusson et al., “Parametric estimation of 3D tubular structures for diffuse optical tomography,” Biomed. Opt. Express, 4 (2), 271 –286 (2013). https://doi.org/10.1364/BOE.4.000271 BOEICL 2156-7085 Google Scholar

35.

L. Wu et al., “Shape-parameterized diffuse optical tomography holds promise for sensitivity enhancement of fluorescence molecular tomography,” Biomed. Opt. Express, 5 (10), 3640 –3659 (2014). https://doi.org/10.1364/BOE.5.003640 BOEICL 2156-7085 Google Scholar

36.

M. E. Kilmer et al., “A shape-based reconstruction technique for DPDW data,” Opt. Express, 7 (13), 481 –491 (2000). https://doi.org/10.1364/OE.7.000481 OPEXFF 1094-4087 Google Scholar

37.

A. Aghasi, M. Kilmer and E. L. Miller, “Parametric level set methods for inverse problems,” SIAM J. Imaging Sci., 4 (2), 618 –650 (2011). https://doi.org/10.1137/100800208 Google Scholar

38.

N. Naik et al., “Fully nonlinear SP₃ approximation based fluorescence optical tomography,” IEEE Trans. Med. Imaging, 36 (11), 2308 –2318 (2017). https://doi.org/10.1109/TMI.2017.2718028 ITMID4 0278-0062 Google Scholar

39.

N. Naik, R. Beatson and J. Eriksson, “Radial-basis-function level-set-based regularized Gauss–Newton-filter reconstruction scheme for dynamic shape tomography,” Appl. Opt., 53 (29), 6872 –6884 (2014). https://doi.org/10.1364/AO.53.006872 APOPAI 0003-6935 Google Scholar

40.

X. Zhang et al., “Shape-based reconstruction of dynamic fluorescent yield with a level set method,” Biomed. Eng. Online, 15 (1), 6 (2016). https://doi.org/10.1186/s12938-016-0124-y Google Scholar

41.

O. Gottam, N. Naik and S. Gambhir, “Shape based pharmacokinetic fluorescence optical tomography,” Proc. SPIE, 10711 107110M (2018). https://doi.org/10.1117/12.2319358 PSISDG 0277-786X Google Scholar

42.

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

43.

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 (16), 2940 –2950 (2003). https://doi.org/10.1364/AO.42.002940 APOPAI 0003-6935 Google Scholar

44.

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 (10), 1861 –1871 (2006). https://doi.org/10.1109/TBME.2006.881796 IEBEAX 0018-9294 Google Scholar

45.

J. A. Jacquez, Compartmental Analysis in Biology and Medicine, Elsevier, Amsterdam (1972). Google Scholar

46.

V. Ntziachristos et al., “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U. S. A., 97 (6), 2767 –2772 (2000). https://doi.org/10.1073/pnas.040570597 Google Scholar

47.

D. A. Hansen et al., “Indocyanine green (ICG) staining and demarcation of tumor margins in a rat glioma model,” Surg. Neurol., 40 (6), 451 –456 (1993). https://doi.org/10.1016/0090-3019(93)90046-4 SGNRAI 0090-3019 Google Scholar

48.

T. Ishizawa et al., “Real-time identification of liver cancers by using indocyanine green fluorescent imaging,” Cancer, 115 (11), 2491 –2504 (2009). https://doi.org/10.1002/cncr.v115:11 CANCAR 0008-543X Google Scholar

49.

M.-Y. Su et al., “Tumor characterization with dynamic contrast–enhanced MRI using mr contrast agents of various molecular weights,” Magn. Reson. Med., 39 (2), 259 –269 (1998). https://doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

50.

S. S. Jambhekar and P. J. Breen, Basic Pharmacokinetics, 76 Pharmaceutical Press, London, United Kingdom (2009). Google Scholar

51.

D. L. Buckley, “Uncertainty in the analysis of tracer kinetics using dynamic contrast-enhanced T1-weighted MRI,” Magn. Reson. Med., 47 (3), 601 –606 (2002). https://doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

52.

P. S. Tofts et al., “Estimating kinetic parameters from dynamic contrast-enhanced T1-weighted MRI of a diffusable tracer: standardized quantities and symbols,” J. Magn. Reson. Imaging, 10 (3), 223 –232 (1999). https://doi.org/10.1002/(ISSN)1522-2586 Google Scholar

53.

C.-T. Chen, Linear System Theory and Design, Oxford University Press, Inc., New York (1999). Google Scholar

54.

F. Fedele, J. Laible and M. Eppstein, “Coupled complex adjoint sensitivities for frequency-domain fluorescence tomography: theory and vectorized implementation,” J. Comput. Phys., 187 (2), 597 –619 (2003). https://doi.org/10.1016/S0021-9991(03)00150-5 JCTPAH 0021-9991 Google Scholar

55.

H. Wendland, Scattered Data Approximation, 17 Cambridge University Press, Cambridge (2004). Google Scholar

56.

T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Trans. Image Process., 10 (2), 266 –277 (2001). https://doi.org/10.1109/83.902291 IIPRE4 1057-7149 Google Scholar

57.

D. Simon, Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches, John Wiley & Sons, Hoboken, New Jersey (2006). Google Scholar

58.

N. Morrison, Tracking Filter Engineering: The Gauss-Newton and Polynomial Filters, Institution of Engineering and Technology, London, United Kingdom (2013). Google Scholar

59.

M. Schweiger, S. R. Arridge and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol., 50 (10), 2365 –2386 (2005). https://doi.org/10.1088/0031-9155/50/10/013 PHMBA7 0031-9155 Google Scholar

60.

J. Nocedal and S. Wright, Numerical Optimization, Series in Operations Research and Financial Engineering, Springer, New York (2006). Google Scholar

61.

E. de Sturler and M. E. Kilmer, “A regularized Gauss–Newton trust region approach to imaging in diffuse optical tomography,” SIAM J. Sci. Comput., 33 (5), 3057 –3086 (2011). https://doi.org/10.1137/100798181 SJOCE3 1064-8275 Google Scholar

62.

A. R. Conn, N. I. Gould and P. L. Toint, Trust Region Methods, 1 SIAM, Philadelphia (2000). https://doi.org/10.1137/1.9780898719857 Google Scholar

63.

M. J. Eppstein et al., “Biomedical optical tomography using dynamic parameterization and Bayesian conditioning on photon migration measurements,” Appl. Opt., 38 (10), 2138 –2150 (1999). https://doi.org/10.1364/AO.38.002138 APOPAI 0003-6935 Google Scholar

64.

M. Schweiger et al., “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys., 22 (11), 1779 –1792 (1995). https://doi.org/10.1118/1.597634 MPHYA6 0094-2405 Google Scholar

65.

M. J. Eppstein et al., “A comparison of exact and approximate adjoint sensitivities in fluorescence tomography,” IEEE Trans. Med. Imaging, 22 (10), 1215 –1223 (2003). https://doi.org/10.1109/TMI.2003.818165 ITMID4 0278-0062 Google Scholar

66.

F. Larusson, S. Fantini and E. L. Miller, “Parametric level set reconstruction methods for hyperspectral diffuse optical tomography,” Biomed. Opt. Express, 3 (5), 1006 –1024 (2012). https://doi.org/10.1364/BOE.3.001006 BOEICL 2156-7085 Google Scholar

67.

H. K. Kim and A. H. Hielscher, “A PDE-constrained SQP algorithm for optical tomography based on the frequency-domain equation of radiative transfer,” Inverse Probl., 25 (1), 015010 (2008). https://doi.org/10.1088/0266-5611/25/1/015010 INPEEY 0266-5611 Google Scholar

68.

A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys., 220 (1), 441 –470 (2006). https://doi.org/10.1016/j.jcp.2006.07.007 JCTPAH 0021-9991 Google Scholar

Biography

Omprakash Gottam received his MTech degree in electrical engineering from Indian Institute of Technology Kanpur (IIT-Kanpur), India, in 2013. Currently, he is a PhD student in the Department of Electrical Engineering at the IIT-Kanpur. His current research interests include reconstruction algorithms for fluorescence and photoacoustic tomography.

Naren Naik is a member of the faculty at the Department of Electrical Engineering as well as the Center for Lasers and Photonics, at the IIT-Kanpur, India. He is a specialist in tomographic reconstruction and analysis, with an emphasis on limited data and subsurface settings, and has developed frameworks for shape-based and dynamic tomographies. He and his research group work on computational aspects of modeling, reconstruction, analysis and calibration, in optical, fluorescence-optical, photoacoustic, electrical impedance and microwave tomographies in applications of early cancer detection and metabolic imaging, subsurface imaging, remote-sensing and battlefield surveillance.

Sanjay Gambhir is head of the Nuclear Medicine Department at Sanjay Gandhi Post Graduate Institute of Medical Sciences (SGPGIMS), Lucknow, India. SGPGIMS is one of the large university teaching hospitals in the state of Uttar-Pradesh. He runs a teaching and research program in the use of radio-nuclides in imaging along with a cyclotron and radionuclide therapy unit. He has research interests in exploring new radio-ligands and combining various radio-nuclidic and nonradionuclidic imaging techniques for clinical applications.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Omprakash Gottam, Naren Naik, and Sanjay Gambhir "Parameterized level-set based pharmacokinetic fluorescence optical tomography using the regularized Gauss–Newton filter," Journal of Biomedical Optics 24(3), 031010 (10 October 2018). https://doi.org/10.1117/1.JBO.24.3.031010

Received: 14 June 2018; Accepted: 14 September 2018; Published: 10 October 2018

Access the abstract

JOURNAL ARTICLE
17 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY