## 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.^{1}2.3.4.5.6.7.8.9.10.11.12.13.^{–}^{14}

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.^{15}16.^{–}^{17} These rates have been found to bear the signatures of abnormalities in the tissue,^{1}^{,}^{4}^{,}^{9}^{,}^{10}^{,}^{18}19.20.21.^{–}^{22} 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,^{9}^{,}^{10}^{,}^{19}20.21.^{–}^{22} 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^{9} and Wang et al.^{10} 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.^{9} comparing the performance of “one-step” linear as well as nonlinear inversions, with the “two-step”^{21} 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.^{13}^{,}^{23}24.^{–}^{25} The structural priors better specify tissue optical properties corresponding to the organs housing various regions in the image, as well as regularizing the reconstructions.^{26} In such a dual-mode setting,^{12} 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^{27}28.29.30.31.32.33.34.^{–}^{35} 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.^{36} 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.^{28} 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.^{29} In an implicit parametrized level-set framework, Aghasi et al.^{37} 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.^{38} 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^{39} 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^{40} 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.^{40} 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.^{41} 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^{41} 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^{9}^{,}^{19}^{,}^{21}^{,}^{42}43.^{–}^{44} 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.^{16}^{,}^{45} Indocyanine green (ICG) is an optical contrast agent, widely used for cancer detection studies.^{46}47.^{–}^{48} A two compartment model (schematic in Fig. 1) has been reported to be suitable for describing ICG pharmacokinetics.^{44} ICG administered intravenously into the blood stream binds to plasma proteins (albumin) and acts as a macromolecular agent.^{43}^{,}^{46} Due to high leaky vasculature in the tumor region,^{43}^{,}^{49} 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.

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.^{50} 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.^{43}^{,}^{44}^{,}^{51} In the tumor regions, ${v}_{e}$ is in the range of 0.2 to 0.5,^{52} whereas ${v}_{p}$ is in the range of 0.013 to 0.067.^{43}^{,}^{44}

Let ${C}_{p}(\mu \mathrm{M})$ and ${C}_{e}(\mu \mathrm{M})$ be the concentration of ICG in the plasma compartment and EES compartment, respectively, ${k}_{pe}({\mathrm{s}}^{-1})$ [respectively, ${k}_{pe}({\mathrm{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}_{elm}({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:^{9}

## (1)

$$\dot{\mathbf{C}}(\overrightarrow{r},t)=\mathbf{K}({k}_{pe}(\overrightarrow{r}),{k}_{ep}(\overrightarrow{r}),{k}_{elm}(\overrightarrow{r}))\mathbf{C}(\overrightarrow{r},t),$$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^{9}^{,}^{53}

## (2)

$$\mathbf{C}(\overrightarrow{r},j+1)=\mathbf{T}({\tau}_{11}(\overrightarrow{r}),{\tau}_{12}(\overrightarrow{r}),{\tau}_{21}(\overrightarrow{r}),{\tau}_{22}(\overrightarrow{r}))\mathbf{C}(\overrightarrow{r},j),$$## (3)

$$\mathbf{T}\equiv {e}^{\mathbf{K}{t}_{s}}\equiv \left[\begin{array}{cc}{\tau}_{11}(\overrightarrow{r})& {\tau}_{12}(\overrightarrow{r})\\ {\tau}_{21}(\overrightarrow{r})& {\tau}_{22}(\overrightarrow{r})\end{array}\right],$$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^{54}

## (4)

$$\begin{array}{l}-\nabla \xb7({D}_{x}\nabla {\mathrm{\varphi}}_{x})+{k}_{x}{\mathrm{\varphi}}_{x}={S}_{x}\\ -\nabla \xb7({D}_{m}\nabla {\mathrm{\varphi}}_{m})+{k}_{m}{\mathrm{\varphi}}_{m}=\beta {\mathrm{\varphi}}_{x}\end{array}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{in}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\Omega},$$## (5)

$$\begin{array}{l}\overrightarrow{n}\xb7({D}_{x}\nabla {\mathrm{\varphi}}_{x})+{b}_{x}{\mathrm{\varphi}}_{x}=0\\ \overrightarrow{n}\xb7({D}_{m}\nabla {\mathrm{\varphi}}_{m})+{b}_{m}{\mathrm{\varphi}}_{m}=0\end{array}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{on}\text{\hspace{0.17em}\hspace{0.17em}}\partial \mathrm{\Omega},$$## (6)

$${D}_{x/m}=\frac{1}{3({\mu}_{a(x/m)i}+{\mu}_{a(x/m)f}+{\mu}_{s(x/m)}^{\prime})},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\beta =\frac{{\varphi}_{q}{\mu}_{axf}}{1-i\omega {\tau}_{l}},\phantom{\rule{0ex}{0ex}}{b}_{x/m}=\frac{1-{R}_{(x/m)}}{2(1+{R}_{(x/m)})},\phantom{\rule[-0.0ex]{1em}{0.0ex}}{k}_{x/m}=\frac{i\omega}{c}+{\mu}_{a(x/m)i}+{\mu}_{a(x/m)f},$$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 $({\mu}_{axf})$ from the measured complex log-intensity at the emission wavelength.^{54} We assume an *a priori* known linear relationship between ${\mu}_{axf}$ and ${\mu}_{amf}$.^{54} The measured log-intensity is as follows:

## (7)

$$\mathrm{log}({\mathrm{\varphi}}_{m}[\{{\overrightarrow{r}}_{d}\}])\equiv \mathrm{log}|{\mathrm{\varphi}}_{m}({\overrightarrow{r}}_{d})|+i\nu ({\overrightarrow{r}}_{d}),$$## (8)

$$y({\overrightarrow{r}}_{d},t)\equiv \mathrm{log}\text{\hspace{0.17em}}{\mathrm{\varphi}}_{m}[\{{\overrightarrow{r}}_{d}\},t]\equiv \mathcal{G}({\underset{\_}{\mu}}_{axf}[\overrightarrow{r},t]),$$^{54}

The relation between the total fluorophore concentration in the region of interest and the absorption coefficient is given by^{9}

## (9)

$${\mu}_{a(x/m)f}(\overrightarrow{r},t)=\mathrm{ln}\text{\hspace{0.17em}}10\xb7{\u03f5}_{(x/m)}\xb7C(\overrightarrow{r},t),$$^{9}

## 2.2.

### Level-Set Representation and State Variable Model

We consider a level-set based representation of a pharmacokinetic parameter ${k}_{\xi}$, with $\xi \in \{ep,pe,elm\}$, as follows:

## (11)

$${k}_{\xi}(\overrightarrow{r})={k}_{\xi}^{i}(\overrightarrow{r})H[{\varphi}_{\gamma}(\overrightarrow{r})]+{k}_{\xi}^{o}(\overrightarrow{r})[1-H({\varphi}_{\gamma}(\overrightarrow{r}))],$$^{37}

The compactly supported RBF level-set function can be written as^{37} follows:

## (12)

$${\varphi}_{\gamma}(\overrightarrow{r})\equiv \varphi (\overrightarrow{r},\underset{\gamma}{\underset{\u23df}{[\mathit{\alpha},\mathit{\zeta},\mathit{\chi}]}})=\sum _{l=1}^{m}{\alpha}_{l}\psi ({\Vert {\zeta}_{l}(\overrightarrow{r}-{\chi}_{l})\Vert}^{\u2020}),$$## (13)

$$\tilde{r}\equiv {\Vert \overrightarrow{r}\Vert}^{\u2020}=\sqrt{{\Vert \overrightarrow{r}\Vert}^{2}+{v}^{2}},$$^{55}$m$ is the number of RBFs used, ${\alpha}_{l}$ is the weighting factor, ${\zeta}_{l}$ is the dilation factor, and ${\chi}_{l}$ denotes the RBF center coordinates.

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

^{55}defined as follows: where $H(x)$ is the Heaviside function. We also use a mollified Heaviside function in our work:

^{56}

## (16)

$${H}_{\u03f5}(\varphi )=\{\begin{array}{ll}1& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}\varphi >{\u03f5}_{w},\\ 0& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}\varphi <-{\u03f5}_{w},\\ \frac{1}{2}+\frac{\varphi}{2{\u03f5}_{w}}+\frac{1}{2\pi}\text{\hspace{0.17em}}\mathrm{sin}\left(\frac{\pi \varphi}{{\u03f5}_{w}}\right)& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}|\varphi |\le {\u03f5}_{w},\end{array}$$## (17)

$${C}_{(e/p)}(\overrightarrow{r},j)={C}_{(e/p)}^{i}(j){H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})]+{C}_{(e/p)}^{o}(j)[1-{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{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:

## (18)

$${v}_{(e/p)}(\overrightarrow{r})={v}_{(e/p)}^{i}(\overrightarrow{r}){H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})]+{v}_{(e/p)}^{o}(\overrightarrow{r})[1-{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r}))].$$Substituting the level-set representation of concentrations and volume fractions in the above equation (for relation between ${\mu}_{axf}$ and $C$), we obtain the following:

## (19)

$$\mu (\overrightarrow{r},j)=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7[({C}_{e}(\overrightarrow{r},j){v}_{e}(\overrightarrow{r})+{C}_{p}(\overrightarrow{r},j){v}_{p}(\overrightarrow{r}))]=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7[({C}_{e}{(j)}^{i}{v}_{e}^{i}+{C}_{p}^{i}(j){v}_{p}^{i}){H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r}))+({C}_{e}^{o}(j){v}_{e}^{o}+{C}_{p}^{o}(j){v}_{p}^{o})(1-{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r})))],$$Using the level-set representation of pharmacokinetic rates and concentrations, the coupled ODE [Eq. (1)] is rewritten as follows:

## (20)

$$\left[\begin{array}{c}{\dot{C}}_{e}^{i}(t)\\ {\dot{C}}_{p}^{i}(t)\end{array}\right]{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r}))+\left[\begin{array}{c}{\dot{C}}_{e}^{o}(t)\\ {\dot{C}}_{p}^{o}(t)\end{array}\right][1-{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r}))]=\left[\begin{array}{cc}-{k}_{ep}^{i}& {k}_{pe}^{i}\\ {k}_{ep}^{i}& -({k}_{pe}^{i}+{k}_{elm}^{i})\end{array}\right]\left[\begin{array}{c}{C}_{e}^{i}(t)\\ {C}_{p}^{i}(t)\end{array}\right]{H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})]+\left[\begin{array}{cc}-{k}_{ep}^{o}& {k}_{pe}^{o}\\ {k}_{ep}^{o}& -({k}_{pe}^{o}+{k}_{elm}^{o})\end{array}\right]\left[\begin{array}{c}{C}_{e}^{o}(t)\\ {C}_{p}^{o}(t)\end{array}\right][1-{H}_{\u03f5}({\varphi}_{\gamma}(\overrightarrow{r}))].$$The time-discretized version for the above equation can thus be written as follows:

## (21)

$$\left[\begin{array}{c}{C}_{e}^{i}(j+1)\\ {C}_{p}^{i}(j+1)\\ {C}_{e}^{o}(j+1)\\ {C}_{p}^{o}(j+1)\end{array}\right]=\left[\begin{array}{cccc}{\tau}_{11}^{i}& {\tau}_{12}^{i}& 0& 0\\ {\tau}_{21}^{i}& {\tau}_{22}^{i}& 0& 0\\ 0& 0& {\tau}_{11}^{o}& {\tau}_{12}^{o}\\ 0& 0& {\tau}_{21}^{o}& {\tau}_{22}^{o}\end{array}\right]\left[\begin{array}{c}{C}_{e}^{i}(j)\\ {C}_{p}^{i}(j)\\ {C}_{e}^{o}(j)\\ {C}_{p}^{o}(j)\end{array}\right],$$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,^{9}^{,}^{10}^{,}^{22}^{,}^{57} or deterministic schemes, such as the GN filter.^{39}^{,}^{58} 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:

## (22)

$${\mathrm{\Theta}}_{j}\equiv \{\underset{\mathit{C}}{\underset{\u23df}{{C}_{e}^{i}(j),{C}_{p}^{i}(j),{C}_{e}^{o}(j),{C}_{p}^{o}(j)}},\underset{\mathit{k}}{\underset{\u23df}{{k}_{pe}^{i},{k}_{ep}^{i},{k}_{elm}^{i},{k}_{pe}^{o},{k}_{ep}^{o},{k}_{elm}^{o}}},\underset{\mathit{v}}{\underset{\u23df}{{v}_{e}^{i},{v}_{e}^{o},{v}_{p}^{i},{v}_{p}^{o}}},\underset{\mathit{\gamma}}{\underset{\u23df}{\alpha ,\beta ,{\chi}_{1},{\chi}_{2}}}\}.$$Assuming our state equation is exact, we would need to explicitly reconstruct only ${\mathrm{\Theta}}_{0}$. We can rewrite the state equation as follows:

## (23)

$$\left[\begin{array}{c}{C}_{e}^{i}(j+1)\\ {C}_{p}^{i}(j+1)\\ {C}_{e}^{o}(j+1)\\ {C}_{p}^{o}(j+1)\\ \mathit{k}\\ \mathit{v}\\ \mathit{\gamma}\end{array}\right]=\left[\begin{array}{ccccccc}{\tau}_{11}^{i}& {\tau}_{12}^{i}& 0& 0& 0& 0& 0\\ {\tau}_{21}^{i}& {\tau}_{22}^{i}& 0& 0& 0& 0& 0\\ 0& 0& {\tau}_{11}^{o}& {\tau}_{12}^{o}& 0& 0& 0\\ 0& 0& {\tau}_{21}^{o}& {\tau}_{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}_{4m}\end{array}\right]\left[\begin{array}{c}{C}_{e}^{i}(j)\\ {C}_{p}^{i}(j)\\ {C}_{e}^{o}(j)\\ {C}_{p}^{o}(j)\\ \mathit{k}\\ \mathit{v}\\ \mathit{\gamma}\end{array}\right],$$## (24)

$${\mathrm{\Theta}}_{j+1}=A({\mathrm{\Theta}}_{j})\xb7{\mathrm{\Theta}}_{j}=f({\mathrm{\Theta}}_{j}),$$## (25)

$$A({\mathrm{\Theta}}_{j})\equiv \left[\begin{array}{ccc}{T}^{\mathrm{i}}& 0& 0\\ 0& {T}^{\mathrm{o}}& 0\\ 0& 0& {I}_{4m+10}\end{array}\right];\phantom{\rule[-0.0ex]{2em}{0.0ex}}{T}^{\mathrm{i}/\mathrm{o}}=\left[\begin{array}{cc}{\tau}_{11}^{\mathrm{i}/\mathrm{o}}(\mathbf{k})& {\tau}_{12}^{\mathrm{i}/\mathrm{o}}(\mathbf{k})\\ {\tau}_{21}^{\mathrm{i}/\mathrm{o}}(\mathbf{k})& {\tau}_{22}^{\mathrm{i}/\mathrm{o}}(\mathbf{k})\end{array}\right].$$The discrete-time measurement equation at time $j$ can be formally written from Eq. (8) as follows:

## 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:^{39}^{,}^{58}

## (27)

$${\widehat{\mathrm{\Theta}}}_{0}=\underset{{\mathrm{\Theta}}_{0}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}\text{\hspace{0.17em}}\mathcal{F}({\mathrm{\Theta}}_{0})\u2254\frac{1}{2}{\Vert (\mathbf{g}({\mathrm{\Theta}}_{0})-\mathbf{y})\Vert}^{2}+\tau R({\mathrm{\Theta}}_{0}),$$$R({\mathrm{\Theta}}_{0})$, the regularization functional, is chosen here as a minimum-norm penalty:

where ${\mathrm{\Theta}}_{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.^{39}

^{,}

^{59}

^{,}

^{60}A regularized GN update ${p}_{\mathrm{\Theta}}$ solves at the current iterate $\mathrm{\Theta}$:

## (29)

$${\widehat{p}}_{\mathrm{\Theta}}=\underset{{p}_{\mathrm{\Theta}}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}{\Vert \begin{array}{c}\mathbf{J}(\mathrm{\Theta}){p}_{\mathrm{\Theta}}+\mathbf{r}\\ \sqrt{\tau}(\mathrm{\Theta}-{\mathrm{\Theta}}_{c}+{p}_{\mathrm{\Theta}})\end{array}\Vert}^{2},$$## (30)

$$\mathbf{J}=\left[\begin{array}{c}{J}_{0}\\ \vdots \\ {J}_{M-1}\end{array}\right],\phantom{\rule[-0.0ex]{2em}{0.0ex}}\mathbf{r}=\left[\begin{array}{c}{r}_{0}\\ \vdots \\ {r}_{M-1}\end{array}\right],$$## (31)

$${J}_{j-1}={G}_{j-1}[{\mathrm{\Theta}}_{j-1}]{F}_{j-1}[{\mathrm{\Theta}}_{j-1}]{F}_{0}[{\mathrm{\Theta}}_{0}],$$## 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,^{32} 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 ${\mathbf{J}}_{\mathbf{a}}{p}_{\mathrm{\Theta}}=-{\mathbf{r}}_{\mathbf{a}}$, where we use the augmented Jacobian ${\mathbf{J}}_{a}\equiv [\begin{array}{c}\mathbf{J}\\ \sqrt{\tau}I\end{array}]$ and the augmented residual ${\mathbf{r}}_{a}\equiv [\begin{array}{c}\mathbf{r}\\ \sqrt{\tau}({\mathrm{\Theta}}_{0}-{\mathrm{\Theta}}_{c})\end{array}]$.

The step ${p}_{\mathrm{\Theta}}$ 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:

## (33)

$$\mathrm{\Delta}={\Vert {({\mathbf{J}}_{a}^{T}{\mathbf{J}}_{a}+\lambda I)}^{-1}{\mathbf{J}}_{a}^{T}{\mathbf{r}}_{a}\Vert}_{2},$$^{60}

^{,}

^{61}

To ensure proper scaling, the variables are scaled as $\mathrm{\Theta}=S\xb7\tilde{\mathrm{\Theta}}$ using a diagonal scaling matrix $S$, with elements given by^{59}

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

## (35)

$$\tilde{\mathrm{\Delta}}={\Vert {({\tilde{\mathbf{J}}}_{a}^{T}{\tilde{\mathbf{J}}}_{a}+\lambda I)}^{-1}{\tilde{\mathbf{J}}}_{a}^{T}{\mathbf{r}}_{a}\Vert}_{2}.$$The update in the scaled domain ${\tilde{p}}_{\mathrm{\Theta}}$ is then defined to satisfy the equation as follows:

## (36)

$$\left[\begin{array}{c}{\tilde{\mathbf{J}}}_{\mathbf{a}}\\ \sqrt{\lambda}I\end{array}\right]{\tilde{p}}_{{\mathrm{\Theta}}_{0}}=-\left[\begin{array}{c}{\mathbf{r}}_{\mathbf{a}}\\ \mathbf{0}\end{array}\right].$$The update in the original domain ${p}_{\mathrm{\Theta}}$ is thus given by

We then calculate the cost function $\mathcal{F}$ and reduction ratio, $\rho =\frac{\text{actual reduction}}{\text{predicted reduction}}$ at the nominal estimate ${\mathrm{\Theta}}^{t}$. The regularization parameter, $\tau $, is reduced if $\rho $ is above a threshold ${\rho}_{th}$ with a minimum limit of ${\tau}_{\mathrm{min}}$. The trust region radius updating rule, values of ${\eta}_{1}$, ${\eta}_{2}$, and parameter ${\gamma}_{\text{bad}}$ are based on the practical algorithm in Conn et al.^{62} The algorithmic flow of reconstruction is shown in Algorithm 1.

## Algorithm 1

Trust region-based iterative regularized GN filter.

1: Initialization: ${\mathrm{\Theta}}^{0}$, ${\mathrm{\Theta}}_{c}={\mathrm{\Theta}}^{0}$, ${\mathrm{\Delta}}_{0}$, ${\eta}_{1}=0.01$, ${\eta}_{2}=0.9$, ${\tau}_{0}=0.8$, $i=0$ |

2: while$i<{i}_{\text{max}}$do |

3: Calculate ${\mathbf{J}}^{i}$, ${\mathbf{r}}^{i}$, $\mathcal{F}({\mathrm{\Theta}}^{i})$ using ${\mathrm{\Theta}}^{i}$ |

4: Calculate $\lambda $ using Eq. (33) |

5: Solve for ${p}_{\mathrm{\Theta}}$ using Eqs. (36) and (37) |

6: ${\mathrm{\Theta}}^{t}={\mathrm{\Theta}}^{i}+{p}_{\mathrm{\Theta}}$ |

7: Evaluate $\mathcal{F}({\mathrm{\Theta}}^{t})$ and $\rho $ |

8: if$\rho >{\eta}_{1}$then |

9: Accept update. ${\mathrm{\Theta}}^{i+1}={\mathrm{\Theta}}^{t}$; |

10: if$\rho >{\rho}_{th}$then |

11: ${\tau}_{i+1}=\mathrm{max}({\tau}_{i}/3,{\tau}_{\mathrm{min}})$; |

12: end if |

13: else |

14: ${\mathrm{\Theta}}^{i+1}={\mathrm{\Theta}}^{i}$; |

15: end if |

16: if$\rho >{\eta}_{2}$then |

17: ${\mathrm{\Delta}}_{i+1}=\mathrm{max}(2.5\xb7\Vert {\tilde{p}}_{\mathrm{\Theta}}\Vert ,{\mathrm{\Delta}}_{i})$; |

18: else if$\rho \ge {\eta}_{1}\&\rho <{\eta}_{2}$then |

19: ${\mathrm{\Delta}}_{i+1}={\mathrm{\Delta}}_{i}$; |

20: else if$\rho \ge 0\&\rho <{\eta}_{1}$then |

21: ${\mathrm{\Delta}}_{i+1}=0.25\xb7\Vert {\tilde{p}}_{\mathrm{\Theta}}\Vert $; |

22: else if$\rho <0$then |

23: ${\mathrm{\Delta}}_{i+1}=\mathrm{min}(0.25\xb7\Vert {\tilde{p}}_{\mathrm{\Theta}}\Vert ,\mathrm{max}(0.0625,{\gamma}_{bad})\xb7{\mathrm{\Delta}}_{i})$; |

24: end if |

25: $i=i+1$; |

26: end while |

27: Choose the stopping iterate when ${\u03f5}_{rel}<tol$ or the data-residual stays stable or toggles across iterations. |

The following relative measure^{39} 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:

## (38)

$${\u03f5}_{\mathrm{rel}}=\frac{\Vert {\mathcal{P}}_{{\mathbf{J}}_{a}}{\mathbf{r}}_{a}\Vert}{\Vert {\mathbf{r}}_{a}\Vert},$$## 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

## (39)

$${G}_{j}=\left[\begin{array}{cccc}\underset{D\times 4}{\underset{\u23df}{\frac{\partial {y}_{j}}{\partial \mathcal{C}(j)}}}& \underset{D\times 6}{\begin{array}{c}\underset{\u23df}{\frac{\partial {y}_{j}}{\partial k}}\end{array}}& \underset{D\times 4}{\begin{array}{c}\underset{\u23df}{\frac{\partial {y}_{j}}{\partial v}}\end{array}}& \underset{D\times 4m}{\begin{array}{c}\underset{\u23df}{\frac{\partial {y}_{j}}{\partial \gamma}}\end{array}}\end{array}\right].$$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}_{pe}^{i},{k}_{ep}^{i},{k}_{elm}^{i},{k}_{pe}^{o},{k}_{ep}^{o},{k}_{elm}^{o},{v}_{e}^{i},{v}_{e}^{o},{v}_{p}^{i},{v}_{p}^{o},\alpha ,\beta ,{\chi}_{1},{\chi}_{2}\}$ is calculated using the chain rule as follows:

## (40)

$$\frac{\partial {y}_{j}({\overrightarrow{r}}_{d})}{\partial s}=\frac{1}{{\mathrm{\varphi}}_{m}({\overrightarrow{r}}_{d},j)}[\frac{\partial {\mathrm{\varphi}}_{m}({\overrightarrow{r}}_{d},j)}{\partial {\mu}_{\mathrm{axf}}(j)}\times \frac{\partial {\mu}_{\mathrm{axf}}(j)}{\partial s}+\frac{\partial {\mathrm{\varphi}}_{m}({\overrightarrow{r}}_{d},j)}{\partial {\mu}_{\mathrm{amf}}(j)}\times \frac{\partial {\mu}_{\mathrm{amf}}(j)}{\partial s}].$$The derivative of the fluence ${\mathrm{\varphi}}_{m}$ with respect to ${\mu}_{\mathrm{axf}}$ is given in Fedele et al.^{54} To derive the sensitivity of ${\mu}_{\mathrm{axf}}$ and ${\mu}_{\mathrm{amf}}$ with respect to unknowns $\mathrm{\Theta}$, consider Eq. (19) for the time-varying absorption coefficient:

## (41)

$$\mu (\overrightarrow{r},j)=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7[({C}_{e}^{i}(j){v}_{e}^{i}+{C}_{p}^{i}(j){v}_{p}^{i}){H}_{\u03f5}(\varphi (\overrightarrow{r}))+({C}_{e}^{o}(j){v}_{e}^{o}+{C}_{p}^{o}(j){v}_{p}^{o})(1-{H}_{\u03f5}(\varphi (\overrightarrow{r})))].$$The sensitivity of $\mu (\overrightarrow{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:

## (42)

$$\frac{\partial \mu (\overrightarrow{r},j)}{\partial {C}_{e/p}^{i}(j)}=\mathrm{ln}10\xb7\u03f5\xb7{v}_{e/p}^{i}{H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})];\text{\hspace{0.17em}}\frac{\partial \mu (\overrightarrow{r},j)}{\partial {v}_{e/p}^{i}}=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7{C}_{e/p}^{i}(j){H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})].$$## (43)

$$\frac{\partial \mu (\overrightarrow{r},j)}{\partial {C}_{e/p}^{o}(j)}=\mathrm{ln}10\xb7\u03f5\xb7{v}_{e/p}^{o}\{1-{H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})]\};\text{\hspace{0.17em}}\frac{\partial \mu (\overrightarrow{r},j)}{\partial {v}_{e/p}^{o}}=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7{C}_{e/p}^{o}(j)\{1-{H}_{\u03f5}[{\varphi}_{\gamma}(\overrightarrow{r})]\}.$$The variation of $\mu (\overrightarrow{r},j)$ with respect to $\mathit{\gamma}\in \{\mathit{\alpha},\mathit{\zeta},\mathit{\chi}\}$ is given by

## (44)

$$\frac{\partial \mu (\overrightarrow{r},j)}{\partial \gamma}=\frac{\partial \mu (\overrightarrow{r},j)}{\partial {C}_{e}(\overrightarrow{r},j)}\frac{\partial {C}_{e}(\overrightarrow{r},j)}{\partial \gamma}+\frac{\partial \mu (\overrightarrow{r},j)}{\partial {C}_{p}(\overrightarrow{r},j)}\frac{\partial {C}_{p}(\overrightarrow{r},j)}{\partial \gamma}+\frac{\partial \mu (\overrightarrow{r},j)}{\partial {v}_{e}(\overrightarrow{r})}\frac{\partial {v}_{e}(\overrightarrow{r})}{\partial \gamma}+\frac{\partial \mu (\overrightarrow{r},j)}{\partial {v}_{p}(\overrightarrow{r})}\frac{\partial {v}_{p}(\overrightarrow{r})}{\partial \gamma},$$## (45)

$$\frac{\partial \mu (\overrightarrow{r},j)}{\partial {C}_{e/p}(\overrightarrow{r},j)}=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7{v}_{e/p}(\overrightarrow{r}),\frac{\partial \mu (\overrightarrow{r},j)}{\partial {v}_{e/p}(\overrightarrow{r})}=\mathrm{ln}\text{\hspace{0.17em}}10\xb7\u03f5\xb7{C}_{e/p}(\overrightarrow{r},j).$$Further, from the level-set representation of ${C}_{e/p}$ Eq. (17) and ${v}_{e/p}$ Eq. (18), we have

## (46)

$$\frac{\partial {C}_{e/p}(\overrightarrow{r},j)}{\partial \gamma}=[{C}_{e/p}^{i}(j)-{C}_{e/p}^{o}(j)]{H}_{\u03f5}^{\prime}[{\varphi}_{\gamma}(\overrightarrow{r})]\frac{\partial \varphi}{\partial \gamma},\frac{\partial {v}_{e/p}(\overrightarrow{r})}{\partial \gamma}=({v}_{e/p}^{i}-{v}_{e/p}^{o}){H}_{\u03f5}^{\prime}[{\varphi}_{\gamma}(\overrightarrow{r})]\frac{\partial \varphi}{\partial \gamma},$$^{37}

## (47)

$$\frac{\partial \varphi}{\partial {\alpha}_{l}}=\psi ({\Vert {\zeta}_{j}(\mathbf{x}-{\chi}_{j})\Vert}^{\u2020}),$$## (48)

$$\frac{\partial \varphi}{\partial {\zeta}_{l}}={\alpha}_{j}{\zeta}_{l}\frac{{\Vert (\mathbf{x}-{\chi}_{l})\Vert}^{2}}{{\Vert {\zeta}_{l}(\mathbf{x}-{\chi}_{l})\Vert}^{\u2020}}{\psi}^{\prime}({\Vert {\zeta}_{j}(\mathbf{x}-{\chi}_{l})\Vert}^{\u2020}),$$## (49)

$$\frac{\partial \varphi}{\partial {\chi}_{l}^{k}}=-{\alpha}_{l}{\zeta}_{l}^{2}\frac{({x}^{k}-{\chi}_{l}^{k})}{{\Vert {\zeta}_{l}(\mathbf{x}-{\chi}_{l})\Vert}^{\u2020}}{\psi}^{\prime}({\Vert {\zeta}_{l}(\mathbf{x}-{\chi}_{l})\Vert}^{\u2020}),$$## 3.3.2.

#### Frechet derivative of state transition function

The state transition equation is given by

## (50)

$$\left[\begin{array}{c}{C}_{e}^{i}(j+1)\\ {C}_{p}^{i}(j+1)\\ {C}_{e}^{o}(j+1)\\ {C}_{p}^{o}(j+1)\\ \mathit{k}(j+1)\\ \mathit{v}(j+1)\\ \mathit{\alpha}(j+1)\\ \mathit{\zeta}(j+1)\\ \mathit{\chi}(j+1)\end{array}\right]=\left[\begin{array}{c}{b}_{1}^{i}(\mathit{k},j)\\ {b}_{2}^{i}(\mathit{k},j)\\ {b}_{1}^{o}(\mathit{k},j)\\ {b}_{2}^{o}(\mathit{k},j)\\ \mathit{k}(j)\\ \mathit{v}(j)\\ \mathit{\alpha}(j)\\ \mathit{\zeta}(j)\\ \mathit{\chi}(j)\end{array}\right],$$## (51)

$${b}_{1}^{i}(\mathit{k},j)\equiv {\tau}_{11}^{i}(\mathit{k}){C}_{e}^{i}(j)+{\tau}_{12}^{i}(\mathit{k}){C}_{p}^{i}(j),\phantom{\rule{0ex}{0ex}}{b}_{2}^{i}(\mathit{k},j)\equiv {\tau}_{21}^{i}(\mathit{k}){C}_{e}^{i}(j)+{\tau}_{22}^{i}(\mathit{k}){C}_{p}^{i}(j),\phantom{\rule{0ex}{0ex}}{b}_{1}^{o}(\mathit{k},j)\equiv {\tau}_{11}^{o}(\mathit{k}){C}_{e}^{o}(j)+{\tau}_{12}^{o}(\mathit{k}){C}_{p}^{o}(j),\phantom{\rule{0ex}{0ex}}{b}_{2}^{o}(\mathit{k},j)\equiv {\tau}_{21}^{o}(\mathit{k}){C}_{e}^{o}(j)+{\tau}_{22}^{o}(\mathit{k}){C}_{p}^{o}(j).$$The variables ${\tau}_{pq}^{i/o}(p,q=\mathrm{1,2})$ are functions of parameters $\{{k}_{ep}^{i/o},{k}_{pe}^{i/o},{k}_{elm}^{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

## 4.

## Numerical Studies

## 4.1.

### Test-Case and Reconstruction Settings

A computational domain of size $4\times 4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$, with the origin as the center, is considered for our numerical test cases.^{63} 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^{65} ${\mu}_{\mathrm{axi}}=0.031\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\mu}_{\mathrm{ami}}=0.00415\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\mu}_{\mathrm{sx}}^{\prime}=10.95\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\mu}_{\mathrm{sm}}^{\prime}=9.29\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, $\tau =0.56\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ns}$, $\varphi =0.016$, ${R}_{x,m}=0.431$, ${\u03f5}_{x}=\mathrm{130,000}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{M}}^{-1}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\u03f5}_{m}=\mathrm{11,000}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{M}}^{-1}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$.

The pharmacokinetic rates mentioned for IDC and AC^{21} are used to obtain synthetic measurement data. We assume $6.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{M}$ concentration of fluorophore is injected^{43} 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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^{59} 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{M}$ and concentration of the fluorophore in EES compartment (${C}_{e}^{i}$ and ${C}_{e}^{o}$) is assumed to be $0\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{M}$. The algorithm is terminated when ${\u03f5}_{\mathrm{rel}}<\mathrm{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).

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}_{pe}^{i}$ | 0.0687 | 0.0896 | 0.0900 | 0.0900 | 0.0453 | 0.0803 | 0.0659 | 0.07 | 0.34 |

${k}_{pe}^{o}$ | 0.0306 | 0.0282 | 0.0282 | 0.0265 | 0.0315 | 0.0225 | 0.0220 | 0.03 | 0.28 |

${k}_{ep}^{i}$ | 0.0496 | 0.0436 | 0.0583 | 0.0470 | 0.0446 | 0.0586 | 0.0484 | 0.01 | 0.18 |

${k}_{ep}^{o}$ | 0.0166 | 0.0142 | 0.0140 | 0.0206 | 0.0185 | 0.0153 | 0.0157 | 0.02 | 0.24 |

${k}_{elm}^{i}$ | 0.00449 | 0.0070 | 0.0047 | 0.0070 | 0.0041 | 0.0028 | 0.0042 | 0.13 | 0.55 |

${k}_{\mathrm{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 |

## a

The 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}_{pe}^{i}$ | 0.0292 | 0.0225 | 0.0199 | 0.0229 | 0.0223 | 0.0286 | 0.0216 | 0.05 | 0.31 |

${k}_{pe}^{o}$ | 0.0114 | 0.0103 | 0.0136 | 0.0133 | 0.0102 | 0.0071 | 0.0075 | 0.05 | 0.37 |

${k}_{ep}^{i}$ | 0.0158 | 0.0118 | 0.0136 | 0.0117 | 0.0100 | 0.0121 | 0.0100 | 0.08 | 0.36 |

${k}_{ep}^{o}$ | 0.0065 | 0.0060 | 0.0066 | 0.0082 | 0.0069 | 0.0067 | 0.0057 | 0.02 | 0.26 |

${k}_{\mathrm{elm}}^{i}$ | 0.0043 | 0.0038 | 0.0025 | 0.0031 | 0.0039 | 0.0032 | 0.0033 | 0.06 | 0.41 |

${k}_{\mathrm{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 |

## a

The 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,^{38} ${E}_{\mathrm{AP}}$, is defined across all the time instants at which measurements are taken as follows:

## (52)

$${E}_{AP}=\left(\right[\frac{\sum _{j=1}^{M}|{\mu}_{rec}^{i}(j){A}_{rec}-{\mu}_{ac}^{i}(j){A}_{ac}|}{\sum _{j=1}^{M}|{\mu}_{ac}^{i}(j){A}_{ac}|}]/M)\times 100\%,$$The area of an object is given by

where ${a}_{\text{element}}$ is the area of an element (a constant in our studies), ${\chi}_{\text{object}}(\xb7)$ 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## (54)

$${\overline{x}}_{\text{object}}=\frac{\sum _{i,j}{x}_{i}{\chi}_{\text{object}}({x}_{i},{y}_{j})}{{A}_{\text{object}}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\overline{y}}_{\text{object}}=\frac{\sum _{i,j}{y}_{j}{\chi}_{\text{object}}({x}_{i},{y}_{j})}{{A}_{\text{object}}}.$$The Euclidean distance between the centroids of reconstructed object and actual object is given by

## (55)

$${E}_{c}=\sqrt{{({\overline{x}}_{rec}-{\overline{x}}_{ac})}^{2}+{({\overline{y}}_{rec}-{\overline{y}}_{ac})}^{2}}.$$The Dice coefficient^{66} 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

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

## (57)

$${E}_{\mathrm{NMSE}}=\frac{{\Vert {\mathit{X}}^{\mathit{r}}-{\mathit{X}}^{\mathit{a}}\Vert}^{2}}{{\Vert {\mathit{X}}^{\mathit{a}}\Vert}^{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 $\mathit{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.;^{9} they solve the pointwise problem with an EKF), we evaluate NMSE values [in the form of $20\text{\hspace{0.17em}}\mathrm{log}(\mathrm{NMSE})\mathrm{dB}$] for the mapped-pointwise reconstructed images of the pharmacokinetic parameters and volume fractions. Our obtained pointwise NMSEs for ${k}_{pe}$ (${k}_{\text{in}}$ in Alacam et al.^{9}) range from $-31$ to $-16.9\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ and those for ${k}_{ep}$ (${k}_{\text{out}}$ in Alacam et al.^{9}) range from $-31.3$ to $-20.94\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ across data SNR levels. The work of Alacam et al.^{9} reports NMSEs of $-19.77$ and $-18.49\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ for ${k}_{\text{in}}$ and ${k}_{\text{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.^{9}^{,}^{10}

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

^{67}or approximations, such as the simplified spherical harmonics ($S{P}_{n}$) ones.^{38}^{,}^{68}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 ${\tau}_{ij}i,j=\{\mathrm{1,2}\}$ in terms of ${k}_{\xi}$, with $\xi =\{ep,pe,elm\}$, are given by Eqs. (58)–(60) (for ease of notations, we have omitted the superscript $i/o$ on $\tau $ and $k$-related quantities):

## (58)

$${\tau}_{12}=\frac{4{k}_{ep}{k}_{pe}({\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}-{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}})}{4{k}_{ep}\sqrt{{\mathrm{\Xi}}_{1}}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\tau}_{21}=-\frac{{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}-{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}}{\sqrt{{\mathrm{\Xi}}_{1}}},$$## (59)

$${\tau}_{11}=\frac{{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}\sqrt{{\mathrm{\Xi}}_{1}}+{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}\sqrt{{\mathrm{\Xi}}_{1}}-{k}_{elm}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}+{k}_{elm}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}+{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}^{u}}-{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}-{k}_{pe}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}+{k}_{pe}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}}{2\sqrt{{\mathrm{\Xi}}_{1}}},$$## (60)

$${\tau}_{22}=\frac{{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}\sqrt{{\mathrm{\Xi}}_{1}}+{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}\sqrt{{\mathrm{\Xi}}_{1}}+{k}_{elm}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}-{k}_{elm}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}-{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}^{u}}+{k}_{ep}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}+{k}_{pe}{\mathrm{e}}^{{\mathrm{\Gamma}}_{1}}-{k}_{pe}{\mathrm{e}}^{{\mathrm{\Gamma}}_{2}}}{2\sqrt{{\mathrm{\Xi}}_{1}}},$$## (61)

$${\mathrm{\Gamma}}_{1}={\mathrm{\Xi}}_{2}-\frac{{t}_{s}\sqrt{{\mathrm{\Xi}}_{1}}}{2};\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\mathrm{\Gamma}}_{2}=\frac{{t}_{s}\sqrt{{\mathrm{\Xi}}_{1}}}{2}+{\mathrm{\Xi}}_{2},$$## (62)

$${\mathrm{\Xi}}_{1}={k}_{elm}^{2}-2{k}_{elm}{k}_{ep}+2{k}_{elm}{k}_{pe}+{k}_{ep}^{2}+2{k}_{ep}{k}_{pe}+{k}_{pe}^{2},$$## (63)

$${\mathrm{\Xi}}_{2}=-\frac{{k}_{elm}{t}_{s}}{2}-\frac{{k}_{ep}{t}_{s}}{2}-\frac{{k}_{pe}{t}_{s}}{2}.$$The derivatives of ${\tau}_{ij}$ with respect to ${k}_{\xi}$ 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

_{3}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

## 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.