Translator Disclaimer
10 October 2018 Parameterized level-set based pharmacokinetic fluorescence optical tomography using the regularized Gauss–Newton filter
Author Affiliations +
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.12.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.1516.17 These rates have been found to bear the signatures of abnormalities in the tissue,1,4,9,10,1819.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,1920.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 Yazici9 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,2324.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 importance2728.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 solved39 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 yield40 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 work41 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 studies9,19,21,4243.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.4647.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.

Fig. 1

Block diagram of two compartment model.

JBO_24_3_031010_f001.png

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 (ve/p) as the ratio of its volume of distribution (Ve/p) and the total volume of distribution (V=Ve+Vp), we see that the volume fraction of ICG in both compartments is greater in the tumor region.43,44,51 In the tumor regions, ve is in the range of 0.2 to 0.5,52 whereas vp is in the range of 0.013 to 0.067.43,44

Let Cp(μM) and Ce(μM) be the concentration of ICG in the plasma compartment and EES compartment, respectively, kpe(s1) [respectively, kpe(s1)] 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 kelm(s1) 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)

C˙(r,t)=K(kpe(r),kep(r),kelm(r))C(r,t),
where ‘r’ denotes spatial coordinates and C(r,t)=[Ce(r,t)Cp(r,t)]; K(kpe(r),kep(r),kelm(r))=[kep(r)kpe(r)kep(r)(kpe(r)+kelm(r))].

The corresponding discrete time-state model corresponding to time instants tj and tj+1 (indexed as j and j+1, respectively) for Eq. (1) is given by9,53

(2)

C(r,j+1)=T(τ11(r),τ12(r),τ21(r),τ22(r))C(r,j),
where

(3)

TeKts[τ11(r)τ12(r)τ21(r)τ22(r)],
and ts=tj+1tj 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 by54

(4)

·(Dxϕx)+kxϕx=Sx·(Dmϕm)+kmϕm=βϕxin  Ω,
being subject to the Robin boundary conditions:

(5)

n·(Dxϕx)+bxϕx=0n·(Dmϕm)+bmϕm=0on  Ω,
where

(6)

Dx/m=13(μa(x/m)i+μa(x/m)f+μs(x/m)),β=ϕqμaxf1iωτl,bx/m=1R(x/m)2(1+R(x/m)),kx/m=iωc+μa(x/m)i+μa(x/m)f,
where “x/m” stands for either “x” (excitation) or “m” (emission), n is the outward normal to the boundary, Sx(W/cm2) is the excitation source with modulation frequency ω(rad/s), ϕx(W/cm2) is the excitation fluence, ϕm(W/cm2) is the emission fluence, Dx/m is the diffusion coefficient at the excitation/emission wavelength, kx/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 cm1) 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=1, bx/m are Robin boundary coefficients, and Rx/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 (μaxf) from the measured complex log-intensity at the emission wavelength.54 We assume an a priori known linear relationship between μaxf and μamf.54 The measured log-intensity is as follows:

(7)

log(ϕm[{rd}])log|ϕm(rd)|+iν(rd),
where ϕm|ϕm(rd)|eiν and rd denotes a detector location. Hence, we can formally express the discrete-domain measurement equation at a time instant t as follows:

(8)

y(rd,t)logϕm[{rd},t]G(μ_axf[r,t]),
where y(rd,t) represents the vector of measurements at time t, μ_axf(t) is the vector of unknown absorption coefficients on the computational grid at time t, and G(·) represents the measurement operator; in our case, G(·) is evaluated using the finite-element method for the solution of the governing fluorescence diffusion model in Eqs. (4) and (5).54

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

(9)

μa(x/m)f(r,t)=ln10·ϵ(x/m)·C(r,t),
where ϵx/m is the fluorophore extinction coefficient and C(r,t) is the total fluorophore concentration in the tissue and is given by9

(10)

C(r,t)=vp(r)Cp(r,t)+ve(r)Ce(r,t).

2.2.

Level-Set Representation and State Variable Model

We consider a level-set based representation of a pharmacokinetic parameter kξ, with ξ{ep,pe,elm}, as follows:

(11)

kξ(r)=kξi(r)H[ϕγ(r)]+kξo(r)[1H(ϕγ(r))],
where kξi(·) and kξo(·) represent kξ(·) 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 ϕγ(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.37

The compactly supported RBF level-set function can be written as37 follows:

(12)

ϕγ(r)ϕ(r,[α,ζ,χ]γ)=l=1mαlψ(ζl(rχl)),

(13)

r˜r=r2+v2,
where v is a small real number, ψ is a compactly supported RBF,55 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 C2 polynomial of degree 5:37,55

(14)

ψ(r˜)=(1r˜)+4(4r˜+1),
where the (·)+ denotes a cutoff function55 defined as follows:

(15)

(x)+=x·H(x),
where H(x) is the Heaviside function. We also use a mollified Heaviside function in our work:56

(16)

Hϵ(ϕ)={1if  ϕ>ϵw,0if  ϕ<ϵw,12+ϕ2ϵw+12πsin(πϕϵw)if  |ϕ|ϵw,
where ϵw is the half-width of the transition region of the Heaviside. Concentrations in different compartments Ce, Cp, which are dependent on pharmacokinetic rates, are similarly assumed to be piecewise constant, and they can be expressed as follows:

(17)

C(e/p)(r,j)=C(e/p)i(j)Hϵ[ϕγ(r)]+C(e/p)o(j)[1Hϵ(ϕγ(r))].

The volume fractions ve and vp being different in healthy and tumor regions can also be similarly expressed via their inside/outside values denoted as vei/o and vpi/o, respectively, as follows:

(18)

v(e/p)(r)=v(e/p)i(r)Hϵ[ϕγ(r)]+v(e/p)o(r)[1Hϵ(ϕγ(r))].

Substituting the level-set representation of concentrations and volume fractions in the above equation (for relation between μaxf and C), we obtain the following:

(19)

μ(r,j)=ln10·ϵ·[(Ce(r,j)ve(r)+Cp(r,j)vp(r))]=ln10·ϵ·[(Ce(j)ivei+Cpi(j)vpi)Hϵ(ϕγ(r))+(Ceo(j)veo+Cpo(j)vpo)(1Hϵ(ϕγ(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:

(20)

[C˙ei(t)C˙pi(t)]Hϵ(ϕγ(r))+[C˙eo(t)C˙po(t)][1Hϵ(ϕγ(r))]=[kepikpeikepi(kpei+kelmi)][Cei(t)Cpi(t)]Hϵ[ϕγ(r)]+[kepokpeokepo(kpeo+kelmo)][Ceo(t)Cpo(t)][1Hϵ(ϕγ(r))].

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

(21)

[Cei(j+1)Cpi(j+1)Ceo(j+1)Cpo(j+1)]=[τ11iτ12i00τ21iτ22i0000τ11oτ12o00τ21oτ22o][Cei(j)Cpi(j)Ceo(j)Cpo(j)],
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,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)

Θj{Cei(j),Cpi(j),Ceo(j),Cpo(j)C,kpei,kepi,kelmi,kpeo,kepo,kelmok,vei,veo,vpi,vpov,α,β,χ1,χ2γ}.

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

(23)

[Cei(j+1)Cpi(j+1)Ceo(j+1)Cpo(j+1)kvγ]=[τ11iτ12i00000τ21iτ22i0000000τ11oτ12o00000τ21oτ22o0000000I60000000I40000000I4m][Cei(j)Cpi(j)Ceo(j)Cpo(j)kvγ],
where ID denotes identity matrix of size D×D(D{6,4,4m}). The above equation can be written in simplified form as follows:

(24)

Θj+1=A(Θj)·Θj=f(Θj),
where f denotes the nonlinear state transition function and matrix A(Θj) is

(25)

A(Θj)[Ti000To000I4m+10];Ti/o=[τ11i/o(k)τ12i/o(k)τ21i/o(k)τ22i/o(k)].

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

(26)

yjgj(Θj)=gj[fj1(f0(Θ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:39,58

(27)

Θ^0=argminΘ0F(Θ0)12(g(Θ0)y)2+τR(Θ0),
where τ is the regularization parameter, R(·) 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:

(28)

R(Θ0)=Θ0Θc2,
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.39,59,60 A regularized GN update pΘ solves at the current iterate Θ:

(29)

p^Θ=argminpΘJ(Θ)pΘ+rτ(ΘΘc+pΘ)2,
where the Jacobian J and the residual r are given by

(30)

J=[J0JM1],r=[r0rM1],
where M denotes the number of time instants and Jacobian at time instant j is given by

(31)

Jj1=Gj1[Θj1]Fj1[Θj1]F0[Θ0],
where G[·] and F[·] are the Frechet derivatives of measurement and state transition functions, respectively. The residual at time instant j is given by

(32)

rj=gj(Θj)yj=gj[fj1(f0(Θ0))]yj.

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 JapΘ=ra, where we use the augmented Jacobian Ja[JτI] and the augmented residual ra[rτ(Θ0Θc)].

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:

(33)

Δ=(JaTJa+λI)1JaTra2,
where we calculate the parameter λ using a suitable root-finding technique.60,61

To ensure proper scaling, the variables are scaled as Θ=S·Θ˜ using a diagonal scaling matrix S, with elements given by59

(34)

Sii=1jJij2+τ.

Thus, from Eq. (33), we see that the scaled-domain counterparts of Ja and Δ, namely, J˜a=Ja·S and Δ˜ (the trust region radius in the scaled domain), respectively, satisfy

(35)

Δ˜=(J˜aTJ˜a+λI)1J˜aTra2.

The update in the scaled domain p˜Θ is then defined to satisfy the equation as follows:

(36)

[J˜aλI]p˜Θ0=[ra0].

The update in the original domain pΘ is thus given by

(37)

pΘ=Sp˜Θ.

We then calculate the cost function F and reduction ratio, ρ=actual reductionpredicted reduction at the nominal estimate Θt. The regularization parameter, τ, is reduced if ρ is above a threshold ρth 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.62 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: whilei<imaxdo
3:  Calculate Ji, ri, 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ρ>η1then
9:    Accept update. Θi+1=Θt;
10:   ifρ>ρththen
11:    τi+1=max(τi/3,τmin);
12:   end if
13:  else
14:   Θi+1=Θi;
15:  end if
16:  ifρ>η2then
17:   Δi+1=max(2.5·p˜Θ,Δi);
18:  else ifρη1&ρ<η2then
19:   Δi+1=Δi;
20:  else ifρ0&ρ<η1then
21:   Δi+1=0.25·p˜Θ;
22:  else ifρ<0then
23:   Δi+1=min(0.25·p˜Θ,max(0.0625,γbad)·Δi);
24:  end if
25:  i=i+1;
26:  end while
27: Choose the stopping iterate when ϵrel<tol or the data-residual stays stable or toggles across iterations.

The following relative measure39 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)

ϵrel=PJarara,
where PJa is the orthogonal projection onto the range space of Ja.

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)

Gj=[yjC(j)D×4yjkD×6yjvD×4yjγD×4m].

The sensitivity equation with respect to s{Cei(j),Cpi(j),Ceo(j),Cpo(j),kpei,kepi,kelmi,kpeo,kepo,kelmo,vei,veo,vpi,vpo,α,β,χ1,χ2} is calculated using the chain rule as follows:

(40)

yj(rd)s=1ϕm(rd,j)[ϕm(rd,j)μaxf(j)×μaxf(j)s+ϕm(rd,j)μamf(j)×μamf(j)s].

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

(41)

μ(r,j)=ln10·ϵ·[(Cei(j)vei+Cpi(j)vpi)Hϵ(ϕ(r))+(Ceo(j)veo+Cpo(j)vpo)(1Hϵ(ϕ(r)))].

The sensitivity of μ(r,j) with respect to {Cei(j),Cpi(j),Ceo(j),Cpo(j),vei,veo,vpi,vpo} is given below:

(42)

μ(r,j)Ce/pi(j)=ln10·ϵ·ve/piHϵ[ϕγ(r)];μ(r,j)ve/pi=ln10·ϵ·Ce/pi(j)Hϵ[ϕγ(r)].

(43)

μ(r,j)Ce/po(j)=ln10·ϵ·ve/po{1Hϵ[ϕγ(r)]};μ(r,j)ve/po=ln10·ϵ·Ce/po(j){1Hϵ[ϕγ(r)]}.

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

(44)

μ(r,j)γ=μ(r,j)Ce(r,j)Ce(r,j)γ+μ(r,j)Cp(r,j)Cp(r,j)γ+μ(r,j)ve(r)ve(r)γ+μ(r,j)vp(r)vp(r)γ,
where [from Eq. (19)]

(45)

μ(r,j)Ce/p(r,j)=ln10·ϵ·ve/p(r),μ(r,j)ve/p(r)=ln10·ϵ·Ce/p(r,j).

Further, from the level-set representation of Ce/p Eq. (17) and ve/p Eq. (18), we have

(46)

Ce/p(r,j)γ=[Ce/pi(j)Ce/po(j)]Hϵ[ϕγ(r)]ϕγ,ve/p(r)γ=(ve/pive/po)Hϵ[ϕγ(r)]ϕγ,
where the sensitivities of the level-set, ϕ with repect to the shape parameters, are given by37

(47)

ϕαl=ψ(ζj(xχj)),

(48)

ϕζl=αjζl(xχl)2ζl(xχl)ψ(ζj(xχl)),

(49)

ϕχlk=αlζl2(xkχlk)ζl(xχl)ψ(ζl(xχl)),
where χlk is the k’th component of center χl.

3.3.2.

Frechet derivative of state transition function

The state transition equation is given by

(50)

[Cei(j+1)Cpi(j+1)Ceo(j+1)Cpo(j+1)k(j+1)v(j+1)α(j+1)ζ(j+1)χ(j+1)]=[b1i(k,j)b2i(k,j)b1o(k,j)b2o(k,j)k(j)v(j)α(j)ζ(j)χ(j)],
where

(51)

b1i(k,j)τ11i(k)Cei(j)+τ12i(k)Cpi(j),b2i(k,j)τ21i(k)Cei(j)+τ22i(k)Cpi(j),b1o(k,j)τ11o(k)Ceo(j)+τ12o(k)Cpo(j),b2o(k,j)τ21o(k)Ceo(j)+τ22o(k)Cpo(j).

The variables τpqi/o(p,q=1,2) are functions of parameters {kepi/o,kpei/o,kelmi/o}, whose expressions are given in the Appendix. The Frechet derivative, Fj, of the state transition function at time instant j is given by

Fj=[τ11iτ21i00b1ikpeib1ikepib1ikelmi0000τ21iτ22i00b2ikpeib2ikepib2ikelmi000000τ11oτ21o000b1okpeob1okepob1okelmo000τ21oτ22o000b2okpeob2okepob2okelmo00000000000I4m+10].

4.

Numerical Studies

4.1.

Test-Case and Reconstruction Settings

A computational domain of size 4×4  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 by65 μaxi=0.031  cm1, μami=0.00415  cm1, μsx=10.95  cm1, μsm=9.29  cm1, τ=0.56  ns, ϕ=0.016, Rx,m=0.431, ϵx=130,000  M1cm1, ϵm=11,000  M1cm1.

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

JBO_24_3_031010_f002.png

The pharmacokinetic rates mentioned for IDC and AC21 are used to obtain synthetic measurement data. We assume 6.5  μM concentration of fluorophore is injected43 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 intensity59 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 1.5  cm and “B” being a single bean shaped phantom with lateral and longitudinal extents being 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-T138.99I-B139.35
I-T233.2I-B233.47
I-T329.77I-B329.75
T Phantom (AC)SNR (dB)B phantom (AC)SNR (dB)
A-T139.02A-B139.29
A-T232.87A-B233.4
A-T329.35A-B329.83

Initial estimates for pharmacokinetic rates are taken in between healthy and tumor values. Initial fluorophore concentration in plasma compartment (Cpi and Cpo) is assumed to be 6.5  μM and concentration of the fluorophore in EES compartment (Cei and Ceo) 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.

JBO_24_3_031010_f003.png

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.

JBO_24_3_031010_f004.png

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 Cein and Cpin plots in the tumor region. (c, e, g) Right column shows Ceout and Cpout 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.

JBO_24_3_031010_f005.png

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 Cein and Cpin plots in the tumor region. (c, e, g) Right column shows Ceout and Cpout 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.

JBO_24_3_031010_f006.png

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 Cein and Cpin plots in the tumor region. (c, e, g) Right column shows Ceout and Cpout 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.

JBO_24_3_031010_f007.png

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 Cein and Cpin plots in the tumor region. (c, e, g) Right column shows Ceout and Cpout 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.

JBO_24_3_031010_f008.png

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
ParameterTrueI-T1I-T2I-T3I-B1I-B2I-B3Row NMSEMaximal normalized error
Cei00.010.010.016×1040.010.018×105a0.01a
Ceo000.002900001×106a0.0029a
Cpi6.56.496.56.56.496.56.458×1060.0073
Cpo6.56.56.56.56.56.466.55×1060.0059
kpei0.06870.08960.09000.09000.04530.08030.06590.070.34
kpeo0.03060.02820.02820.02650.03150.02250.02200.030.28
kepi0.04960.04360.05830.04700.04460.05860.04840.010.18
kepo0.01660.01420.01400.02060.01850.01530.01570.020.24
kelmi0.004490.00700.00470.00700.00410.00280.00420.130.55
kelmo0.004460.00290.00470.00380.00410.00280.00400.050.38
vei0.30.44230.43900.47700.35560.31440.33850.140.58
veo003.1×10600001.5×1012a3.1×106a
vpi0.06000.07000.04980.07000.06690.07000.07000.020.16
vpo0.02000.01900.02000.01600.01980.01670.01750.010.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
ParameterTrueA-T1A-T2A-T3A-B1A-B2A-B3Row NMSEMaximal normalized error
Cei00.010.010.010.010.001708.37×105a0.01a
Ceo00000000a0a
Cpi6.56.56.56.56.56.496.52×1071.3×103
Cpo6.56.56.56.56.496.56.51×10103×105
kpei0.02920.02250.01990.02290.02230.02860.02160.050.31
kpeo0.01140.01030.01360.01330.01020.00710.00750.050.37
kepi0.01580.01180.01360.01170.01000.01210.01000.080.36
kepo0.00650.00600.00660.00820.00690.00670.00570.020.26
kelmi0.00430.00380.00250.00310.00390.00320.00330.060.41
kelmo0.00350.00350.00250.00310.00390.00320.00330.020.28
vei0.20000.07940.06450.08910.13170.09010.09670.30.67
veo00000000a0a
vpi0.04000.04720.02910.03420.04280.04160.03770.020.27
vpo0.02000.01910.01990.01870.01960.01700.01850.010.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 EAP, is defined across all the time instants at which measurements are taken as follows:

(52)

EAP=([j=1M|μreci(j)Arecμaci(j)Aac|j=1M|μaci(j)Aac|]/M)×100%,
where Aac and Arec represent the area of actual and reconstructed object, respectively, μreci(j) [respectively, μaci(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

(53)

Aobject=aelementi,jχobject(xi,yj),
where aelement is the area of an element (a constant in our studies), χobject(·) 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)

x¯object=i,jxiχobject(xi,yj)Aobject,y¯object=i,jyjχobject(xi,yj)Aobject.

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

(55)

Ec=(x¯recx¯ac)2+(y¯recy¯ac)2.

The Dice coefficient66 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

(56)

D(S,H)=2|SH||S|+|H|,
where |SH| 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 Xr with respect to its actual value Xa is defined as follows:

(57)

ENMSE=XrXa2Xa2.

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

PhantomEAP (%)EC (cm)D(S,H)ENMSE(k)
I-T11.29(0.04, 0.04)0.840.058
I-T22(0.14, 0.07)0.800.064
I-T31.24(0.04, 0.08)0.810.059
I-B10.20.010.880.068
I-B28.9×1022×1030.890.034
I-B30.340.070.790.01
A-T10.47(0.02,0.03)0.590.047
A-T20.66(0.12,0.05)0.460.076
A-T30.97(0.21,0.03)0.50.049
A-B10.510.020.910.063
A-B20.590.0120.830.026
A-B30.850.070.840.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 20log(NMSE)dB] for the mapped-pointwise reconstructed images of the pharmacokinetic parameters and volume fractions. Our obtained pointwise NMSEs for kpe (kin in Alacam et al.9) range from 31 to 16.9  dB and those for kep (kout in Alacam et al.9) range from 31.3 to 20.94  dB across data SNR levels. The work of Alacam et al.9 reports NMSEs of 19.77 and 18.49  dB for kin and kout, 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.

PhantomE(kpe) dBE(kep) dBE(ve) dBE(vp) dBE(kelm) dB
I-T13124.2217.0363.8218.55
I-T230.925.0717.3667.1151.5
I-T326.922.0110.2549.5731.11
I-B13028.926.26544.4
I-B222.628.5738.535916.8
I-B320.7625.7420.135140.05
A-T128.2630.1113.76061.49
A-T222.6124.659.3879.2921.59
A-T323.520.9410.3873.335.18
A-B130.4928.332.0186.139.62
A-B216.9231.319.3162.7339.74
A-B318.5526.1120.279.3842.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 RTE67 or approximations, such as the simplified spherical harmonics (SPn) 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 τiji,j={1,2} in terms of kξ, with ξ={ep,pe,elm}, are given by Eqs. (58)–(60) (for ease of notations, we have omitted the superscript i/o on τ and k-related quantities):

(58)

τ12=4kepkpe(eΓ2eΓ1)4kepΞ1,τ21=kepeΓ1kepeΓ2Ξ1,

(59)

τ11=eΓ1Ξ1+eΓ2Ξ1kelmeΓ1+kelmeΓ2+kepeΓ1ukepeΓ2kpeeΓ1+kpeeΓ22Ξ1,

(60)

τ22=eΓ1Ξ1+eΓ2Ξ1+kelmeΓ1kelmeΓ2kepeΓ1u+kepeΓ2+kpeeΓ1kpeeΓ22Ξ1,
where Ξ1, Ξ2, Γ1, and Γ2 are given by the following expressions:

(61)

Γ1=Ξ2tsΞ12;Γ2=tsΞ12+Ξ2,

(62)

Ξ1=kelm22kelmkep+2kelmkpe+kep2+2kepkpe+kpe2,

(63)

Ξ2=kelmts2kepts2kpets2.

The derivatives of τij 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. 

Meikle S. R. 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. 

Reutter B., Gullberg G., Huesman R., “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. 

Huesman R. H. 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. 

Langer O. 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. 

Kamasak M. E. 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. 

Wang G., Qi J., “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. 

Milstein A. B. 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. 

Godavarty A. 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. 

Alacam B., Yazici B., “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. 

Wang X. 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. 

Zhang X. 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. 

Liu X. 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. 

Zhang G. 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. 

Zhang X. 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. 

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

16. 

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

17. 

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

18. 

Sun H. 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. 

Gurfinkel M. 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. 

Milstein A. B., Webb K. J., Bouman C. A., “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. 

Alacam B. 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. 

Wang X. 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. 

Zhang G. 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. 

Zhang G. 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. 

Zhang G. 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. 

Davis S. C. 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. 

Santosa F., “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. 

Kilmer M. E. 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. 

Zacharopoulos A. D. 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. 

Dorn O., Lesselier D., “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. 

Arridge S. R. 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. 

Naik N. 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. 

Dorn O., Lesselier D., “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. 

Larusson F. 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. 

Wu L. 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. 

Kilmer M. E. 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. 

Aghasi A., Kilmer M., Miller E. L., “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. 

Naik N. et al., “Fully nonlinear SP3 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. 

Naik N., Beatson R., Eriksson J., “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. 

Zhang X. 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. 

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

42. 

Intes X. 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. 

Cuccia D. J. 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. 

Alacam B. 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. 

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

46. 

Ntziachristos V. 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. 

Hansen D. A. 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. 

Ishizawa T. 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. 

Su M.-Y. 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. 

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

51. 

Buckley D. L., “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. 

Tofts P. S. 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. 

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

54. 

Fedele F., Laible J., Eppstein M., “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. 

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

56. 

Chan T. F., Vese L. A., “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. 

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

58. 

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

59. 

Schweiger M., Arridge S. R., Nissilä I., “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. 

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

61. 

de Sturler E., Kilmer M. E., “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. 

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

63. 

Eppstein M. J. 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. 

Schweiger M. 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. 

Eppstein M. J. 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. 

Larusson F., Fantini S., Miller E. L., “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. 

Kim H. K., Hielscher A. H., “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. 

Klose A. D., Larsen E. W., “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.

© 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.
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
JOURNAL ARTICLE
17 PAGES


SHARE
Advertisement
Advertisement
Back to Top