Dynamic Cell Imaging in PET with Optimal Transport Regularization

TL;DR

This paper introduces a new PET image reconstruction method using optimal transport regularization that effectively tracks small numbers of cells with high accuracy, even under low radioactivity conditions, by enforcing temporal consistency.

Contribution

The method uniquely combines all detected events for dynamic reconstruction and automatically determines the number of cells without additional complexity.

Findings

Robust performance on simulated data with low radioactivity

Automatic detection of number of tracked cells

Achieves 5.3 mm tracking accuracy for moving cells

Abstract

We propose a novel dynamic image reconstruction method from PET listmode data that could be particularly suited to tracking single or small numbers of cells. In contrast to conventional PET reconstruction our method combines the information from all detected events not only to reconstruct the dynamic evolution of the radionuclide distribution, but also to improve the reconstruction at each single time point by enforcing temporal consistency. This is achieved via optimal transport regularization where in principle, among all possible temporally evolving radionuclide distributions consistent with the PET measurement, the one is chosen with least kinetic motion energy. The reconstruction is found by convex optimization so that there is no dependence on the initialization of the method. We study its behaviour on simulated data of a human PET system and demonstrate its robustness even in settings with very low radioactivity. In contrast to previously reported cell tracking algorithms, our technique is oblivious to the number of tracked cells. Without any additional complexity one or multiple cells can be reconstructed, and the model automatically determines the number of particles. For instance, four radiolabelled cells moving at a velocity of 3.1 mm/s and a PET recorded count rate of 1.1 cps (for each cell) could be simultaneously tracked with a tracking accuracy of 5.3 mm inside a simulated human body.

Equations

$$E=(E_{i,j}{)}_{i=1,\dots ,M,j=1\dots ,N}.$$

$$\rho =({\rho }_t(x){)}_{t\in [0,T],x\in \Omega }.$$

$$\frac{\log 2}{T_{\frac{1}{2}}}{\int }_{{\tau }_i}{\int }_S{\rho }_t(x)\mathrm{d}x\mathrm{d}t$$

$$A_{i,j}^{\mathrm{d}}\rho =\frac{\log 2}{T_{\frac{1}{2}}}{\int }_{{\tau }_i}{\int }_{{\Omega }^2}{H}_j(y)G(y,x)p^{\mathrm{d}}{\rho }_t(x)\mathrm{d}x\mathrm{d}y\mathrm{d}t.$$

$$A_{i,j}^{\mathrm{s}}\rho =\frac{1}{N}\frac{\log 2}{T_{\frac{1}{2}}}{\int }_{{\tau }_i}{\int }_{\Omega }{p}^{\mathrm{s}}{\rho }_t(x)\mathrm{d}x\mathrm{d}t.$$

$$\omega =({\omega }_t(x){)}_{t\in [0,T],x\in \Omega }.$$

$${\partial }_t{\rho }_t(x)+\mathrm{div}{\omega }_t(x)=0.$$

$$\mathcal{S}(\rho ,\omega )=\{\begin{array}{ll}{\int }_0^T{\int }_{\Omega }\frac{\Vert {\omega }_t{\Vert }_2^2}{{\rho }_t}\mathrm{d}x\mathrm{d}t& \text{if ρ≥0 and eqrefeqn:continuity holds,}\\ \infty & \text{else.}\end{array}$$

$$P(E|\rho )=\prod _{i=1}^M\prod _{j=1}^N{\mathcal{P}}_{K_{i,j}}(E_{i,j}).$$

$$-\log P(E|\rho )$$

$$J^E[\rho ,\omega ]=\sum _{i=1}^M\sum _{j=1}^N\left[A_{i,j}\rho -E_{i,j}\log A_{i,j}\rho \right]+\beta \mathcal{S}(\rho ,\omega )$$

$$J^{E,p}[\rho,\omega]=\sum_{i=1}^{M}\sum_{j=1}^{N}A_{i,j}\rho-E_{i,j}\log\big{[}A_{i,j}^{\mathrm{d}}\rho+pA_{i,j}^{\mathrm{s}}\rho\big{]}\\ +\beta\,\mathcal{S}(\rho,\omega)\,.$$

$$\sum _{i=1}^M\sum _{j=1}^N{A}_{i,j}\rho ={\int }_0^T{\int }_{\Omega }r(x){\rho }_t(x)\mathrm{d}x\mathrm{d}t.$$

$$J^{E,p}[\rho,\omega]=\int_{0}^{T}\int_{\Omega}r(x)\,\rho_{t}(x)\,{\mathrm{d}}x\,{\mathrm{d}}t\\ -\sum_{i=1}^{M}\sum_{j=1}^{N}E_{i,j}\log\big{[}A_{i,j}^{\mathrm{d}}\rho+pA_{i,j}^{\mathrm{s}}\rho\big{]}+\beta\,\mathcal{S}(\rho,\omega),$$

$${\rho }_t$$

$$\mathcal{S}(\rho ,\omega )$$

$$A_{i,j}^{\mathrm{d}}\rho$$

$$J^{E,p}[\rho,\omega]\approx m+\sum_{i=1}^{M}\sum_{j=1}^{N}E_{i,j}\big{[}\tfrac{d_{i,j}^{2}}{2\epsilon^{2}}-\log m\big{]}\\ +\beta\,m\,\int_{0}^{T}\|X^{\prime}(t)\|_{2}^{2}\,{\mathrm{d}}t+{\mathrm{const}}\,.$$

$$-\log (A_{i,j}^{\mathrm{d}}\rho +p{A}_{i,j}^{\mathrm{s}}\rho )$$

$$\underset{(\rho ,\omega )}{\inf }F(K(\rho ,\omega ))+H(\rho ,\omega )$$

$$\displaystyle d_{{\mathrm{WFR}},\alpha}^{2}(\delta_{x},\delta_{y})=16\,\alpha^{2}\,\sin^{2}\big{(}\min\big{\{}\tfrac{\|x-y\|}{4\alpha},\tfrac{\pi}{4}\big{\}}\big{)}$$

$${\mathrm{err}}^2=\frac{1}{T}{\int }_0^T{d}_{\mathrm{WFR},\alpha }^2(\mathcal{N}({\rho }_t),\mathcal{N}({\rho }_t^{\mathrm{gt}}))\mathrm{d}t$$

$${\mathrm{err}}_{\text{disc}}=\Delta L=2.5\text{mm}.$$

$$\underset{x}{\inf }F(Kx)+H(x)$$

$$K$$

$$F({\rho }_{\text{c}},{\omega }_{\text{c}},\tilde{\rho },\mu )$$

$$F_1({\rho }_{\text{c}},{\omega }_{\text{c}})$$

$$F_2(\tilde{\rho })$$

$$F_3(\mu )$$

$$H(\rho ,\omega )$$

$$y^{(\ell +1)}$$

$$x^{(\ell +1)}$$

Full text

I Introduction

Immune and cell-based therapies have become of emerging interest for eradicating or controlling intracellular pathogens that are difficult to treat with conventional therapies [1]. These immunotherapies generally utilize the primary function of the immune system to kill bacterial and viral pathogens or fight cancer cells. Molecular imaging has great potential to detect and identify biomarkers allowing to treat only those patients which most likely respond to a dedicated immunotherapy and, most importantly, to monitor the treatment’s success [2].

Treatment planning and control need a quantitative measure of the biomarker’s distribution within the patient’s body. Positron emission tomography (PET) is a promising molecular imaging technique as it provides quantitative whole-body data at highest molecular sensitivity. The advent of new radiopharmaceuticals based on F-18 or long-lived Zr-89 and Cu-64 (e.g. [3, 4, 5]) led to the development of new immuno-PET imaging approaches visualizing immune cell behaviour in-vivo.

Besides the use of PET for the visualization and quantitation of biomarkers, a new strategy has been developed recently to follow single cells over time using the excellent sensitivity of PET systems in combination with mathematical modelling techniques (see for instance [6] and references therein). In these approaches, radio-labelled cells, directly injected or indirectly labelled, are followed over time which is different to standard PET techniques where usually a broad spectrum marker like 18F-FDG distributes slowly within the body and shows sites of its accumulation. Therefore, the PET information by singular labelled cells is directly related to the process of cellular motion rather than just showing the cell distribution after long-term tracer accumulation. This application is of great interest to understand the underlying biological processes in immunotherapy. In addition, for the first time, the possibility of tracking the path of individual cells in the body opens up which is of interest to better understand the role and function of circulating tumor cells in the development of cancer metastases.

Note that conventional spatiotemporal PET-reconstruction from listmode data (e. g. [7, 8]) rests upon a (temporally or spatially) smooth parameterization of the radionuclide density (e. g. via kernel methods or cubic splines) and typically also a decoupling of space and time parameterization. Both these features are incompatible with spatiotemporal radionuclide densities induced by individually moving labelled cells (they could only be utilized if the moving density is spatially smooth and the motion is slow). Lee et al. instead proposed a novel algorithm to reconstruct the trajectory of a moving cell directly from PET list-mode data [6]. The trajectory is modelled as a 3D B-spline function of time, and the mean squared distance between the trajectory and the recorded coincidence events is minimized via non-linear optimization. This approach allows to track single sources in a small animal PET setup with an accuracy below 3 mm provided that its activity (in Bq) exceeds four times its velocity (in mm/s). Using the same algorithm, Ouyang et al. demonstrated experimentally the feasibility of tracking a moving point source in another small animal PET system based on BGO detectors without intrinsic background radiation [9]. Due to the reduced noise level, the PET system could recognize moving sources with an activity-to-velocity ratio $>3.45$ decay/mm, outperforming the LSO-based PET system. Langford et al. were able to track multiple labelled yeast cells by following local maxima in a discretized, smoothed back-projection of the PET signal, albeit at significantly higher activity-to-velocity ratios [10].

Our proposed method combines the information from all detected events to reconstruct the dynamic evolution and simultaneously improve the reconstruction quality by enforcing temporal consistency. Our approach combines the advantages of [6, 9] and [10]. It can track multiple particles, the number of which is automatically determined and does not affect the complexity, at activity-to-velocity-ratios comparable with [6]. For single particles our functional is closely related to that of [6]. We tested the new algorithm on simulated data of a human PET system to evaluate a potential clinical use of PET based cell tracking under realistic conditions.

II Mathematical reconstruction model

II-A Measurements and material distribution

The interior of the PET scanner (the measurement volume) will be denoted by $\Omega\subset{\mathbb{R}}^{3}$. Measurements are performed over a time interval $[0,T)$ which we divide into $M$ subintervals $(\tau_{i})_{i=1}^{M}$ of duration $\Delta T=T/M$ with $\tau_{i}=[(i-1)\,\Delta T,i\,\Delta T)$.

With $j\in\{1,\ldots,N\}$ we enumerate all scanner detector pairs. The line of response (LOR) $l_{j}$ of pair $j$ is a straight line through $\Omega$ that connects the centers of the two detectors. In an idealized setting any $\beta$-decay that led to the activation of pair $j$ is assumed to have happened on $l_{j}$. More realistically, due to a variety of effects, the $\beta$-decay either happened close to the line or was caused by a random incidence or strong scattering. The number of detected events in time interval $\tau_{i}$ and detector pair $j$ is denoted $E_{i,j}$, and we abbreviate

[TABLE]

Our goal will be to reconstruct from the measurement $E$ the corresponding radionuclide density $\rho_{t}(x)$ at each position $x\in\Omega$ and time $t\in[0,T]$. We write this as

[TABLE]

The expected number of decays in time interval $\tau_{i}$ and a spatial domain $S\subset\Omega$ is given by

[TABLE]

where ${T_{\frac{1}{2}}}$ denotes the radionuclide halflife.

II-B Forward operator

Each event count $E_{i,j}$ is Poisson distributed with mean $K_{i,j}$ which can be computed from $\rho$ via the forward operator. In our simplified model we distinguish three different outcomes of a positron decay and subsequent photon emission:

1. a)

Detection: The photons undergo at most minor scattering before being registered by a pair of detectors. 2. b)

Scattering: At least one of the photons undergoes substantial scattering which significantly alters its direction. 3. c)

No detection: The emitted photon pair is not detected, e.g. due to absorption or imperfect scanner sensitivity.

The outcomes of the former two will be produced by linear operators $A^{\mathrm{d}}_{i,j}$ and $A^{\mathrm{s}}_{i,j}$; we set $A_{i,j}=A^{\mathrm{d}}_{i,j}+A^{\mathrm{s}}_{i,j}$ and $K_{i,j}=A_{i,j}\rho$. The latter does not lead to any detected events and thus does not contribute to $K_{i,j}$. Note that for each $(i,j)$, $A_{i,j}$ is a linear function that takes the spatiotemporal material distribution $\rho$ as input and outputs the expected number of detected events in time interval $\tau_{i}$ and detector pair $j$. We now describe the two operators.

Detection

We set

[TABLE]

$H_{j}(y)$ encodes the tomography geometry and gives the probability that a collinear unscattered photon pair emitted at $y$ activates the detector pair $j$. $G(y,x)$ accounts for the positron range and gives the probability that a positron emitted at $x$ annihilates at $y$. $G$ can also be used to approximately model minor scattering or slight non-collinearity of the photon pair. We will just consider a spatially homogeneous Gaussian kernel $G(y,x)=\exp(-|y-x|^{2}/(2\epsilon^{2}))/(\sqrt{2\pi}\epsilon)^{3}$ for some fixed $\epsilon>0$. $p^{\mathrm{d}}\in[0,1]$ denotes the probability that a photon pair is neither attenuated nor scattered substantially.

Without complicating the approach one could make $p^{\mathrm{d}}$ space-dependent and use $H_{j}$ and $G$ to encode spatially varying material properties and scanner sensitivity, however, for simplicity our simulations are performed for constant sensitivity.

Scattering

For simplicity we assume that scattering changes the photon rays randomly such that the probability of arriving at a detector pair $j$ is homogeneous,

[TABLE]

Similar as above, $p^{\mathrm{s}}\in[0,1]$ is the probability for strong scattering. As with $A^{\mathrm{d}}_{i,j}$, one can also consider a more elaborate operator that accounts for spatially inhomogeneous scattering. However, we will later see that for our data this approximation is appropriate (cf. fig. 2 and corresponding text). With probability $1-p^{\mathrm{d}}-p^{\mathrm{s}}>0$ photons remain undetected, e.g. due to attenuation. Random coincidences can be modelled by adding a constant background to $K_{i,j}$.

II-C Optimal transport regularization

Our method aims at a regime of low radiation activity with relatively few detected events. Reconstructing the particle distribution $\rho$ from the measurements $E$ is therefore a classical underdetermined and ill-posed inverse problem which we must regularize by incorporating additional prior knowledge. If the radionuclide distribution $\rho$ were temporally almost constant, then by taking long enough time intervals $\tau_{i}$ one would detect enough decays within $\tau_{i}$ to reconstruct the sought $\rho$. However, for quickly changing radionuclide distributions such a framewise reconstruction separately for each time interval $\tau_{i}$ immediately breaks down. Instead it becomes necessary to combine information from multiple consecutive time intervals, of course taking into account that the radionuclide distribution changes between the time intervals.

As $\rho$ describes the physical motion of radioactive particles there will be some temporal consistency, and not all spatiotemporal distributions $\rho$ are equally likely. In our variational reconstruction approach we model this by a penalty term that prefers physically plausible distributions. For this we consider the velocity $v_{t}(x)\in{\mathbb{R}}^{3}$ of the radioactive material at each position $x$ and time $t$, and we will penalize too high velocities. Mathematically, though, it is easier to work with the mass flux $\omega_{t}(x)=\rho_{t}(x)v_{t}(x)$ instead of the velocity, and we abbreviate

[TABLE]

The mass flux satisfies the classical mass conservation

[TABLE]

Here we assume $T\ll{T_{\frac{1}{2}}}$ so that the amount of radioactive material does essentially not change over time; otherwise the constraint could easily be complemented with an additional decay term (leading to so-called unbalanced transport, cf. [11]). In the end we will not only reconstruct $\rho$, but also $\omega$.

To quantify the likelihood of a path $(\rho,\omega)$ we associate with it its physical action

[TABLE]

Note that $\omega_{t}$ can be interpreted as physical momentum of the motion of $\rho_{t}$, so $\int_{\Omega}\tfrac{\|\omega_{t}\|_{2}^{2}}{\rho_{t}}\,{\mathrm{d}}x$ is essentially the kinetic energy of all particles at time $t$. This action is the so-called Benamou–Brenier functional for optimal transport, and the minimum value that can be achieved for a transport of $\rho_{0}$ to $\rho_{T}$ is the squared Wasserstein-2 distance between both measures [12]. We will below add $\beta\cdot\mathcal{S}(\rho,\omega)$ as a penalty term to our log-likelihood functional where $\beta\geq 0$ is a weighting parameter. In essence, we assign a lower penalty to a path $(\rho,\omega)$ the less the mass moves. This acts as a temporal regularization of the particle trajectories (a spatial regularization could be added on top, but since we do not want to impose any further structure on the moving radiotracer distribution this is not done here). Note that $\mathcal{S}$ is a convex function of $(\rho,\omega)$ and thus amenable to global optimization (see [12] for details).

II-D Penalized log-likelihood functional and unbiasing

Denoting by $\mathcal{P}_{\lambda}(k)=\lambda^{k}\mathrm{e}^{-\lambda}/k!$ the Poisson distribution with mean $\lambda$ and writing as before $K_{i,j}=A_{i,j}\rho$, the conditional probability of having the measurement $E$ if the radioactive mass distribution is $\rho$ can be calculated as

[TABLE]

The corresponding negative log-likelihood is

[TABLE]

where the terms $\log E_{i,j}!$ are subsumed into a constant. Inserting $K_{i,j}=A_{i,j}\rho$ and adding the kinetic regularization leads to a penalized maximum likelihood (ML) functional

[TABLE]

to be minimized for $\rho$ and $\omega$.

Unfortunately, the ML estimate in this setting is strongly biased towards declaring every detected photon pair as being unscattered. Essentially, this is due to $A^{\mathrm{s}}_{i,j}\rho$ being very small compared to $A^{\mathrm{d}}_{i,j}\rho$ in the definition of $A_{i,j}$ (indeed, $A^{\mathrm{s}}_{i,j}$ distributes the intensity evenly over all detectors, while $A^{\mathrm{d}}_{i,j}$ concentrates all the intensity on a few detectors, thus producing much higher intensities there). This is a well-known deficiency of the ML estimator. A simple and intuitive way to compensate for this bias is by introducing a tuning parameter $p\geq 0$ that reweighs the contribution of the two forward operators in that term which contains the detections $E_{i,j}$,

[TABLE]

In the following we will assume that events $E_{i,j}$ were most likely caused by unscattered photons if $A_{i,j}^{\mathrm{d}}\rho>pA_{i,j}^{\mathrm{s}}\rho$ and by scattering otherwise. This interpretation can be justified by a more detailed derivation of the transition from (3) to (4). For the present article the intuition will suffice that the unscattered contribution $A_{i,j}^{\mathrm{d}}\rho$ dominates for $p$ close to zero, whereas the scattering contribution $pA_{i,j}^{\mathrm{s}}\rho$ dominates for very large $p$. By choosing an intermediate $p$, just the right amount of events will be interpreted as scattered.

Finally, recall that we will be working in the regime of very low radiation activity which means that $E$ is expected to be sparse, i.e. only a small fraction of entries $E_{i,j}$ will be non-zero. $\sum_{i=1}^{M}\sum_{j=1}^{N}A_{i,j}\rho$ is the expected total number of detected events (with and without scattering). Consequently, there is a function $r:\Omega\to{\mathbb{R}}$, related to the material properties and the sensitivity of the PET scanner, such that

[TABLE]

One can then rewrite (4) as

[TABLE]

where the sum now only runs over the sparse subset where $E_{i,j}>0$, and only those elements of the full tomography operator must be computed. In our experiments the area where the particles are moving is small compared to the dimensions of the PET scanner. Thus we may assume that $r$ is spatially constant. Moreover, one can show that if $(\rho,\omega)$ minimize $J^{E,p}$ for some $r$ and $\beta$ then $(\rho,\omega)/q$ will minimize $J^{E,p}$ for $q\cdot r$ and $q\cdot\beta$, for any $q>0$. Hence, only the relative scaling of $r$ and $\beta$ is relevant, and we may for simplicity set $r(x)=1/T$.

II-E Relation to single-cell tracking functional

At first glance our functional looks very different from the tracking functional for a single cell by Lee et al. in[6]. In this section we illustrate that there is a close relation when our functional is constrained to a single particle; in fact, the functional by Lee et al. is a special case of our restricted functional up to slight implementation and modelling variations. The advantage of our new formulation is the straightforward and easy applicability to tracking of multiple cells (with potentially different unknown amounts of radioactivity).

Let $t\mapsto X(t)$ describe the trajectory of a single particle in $\Omega$; for simplicity assume that $X$ is differentiable such that $X^{\prime}(t)$ denotes the particle velocity. The particle mass is denoted by $m\geq 0$. Then the corresponding time-dependent particle distribution $\rho$ and the material flux $\omega$ are

[TABLE]

where $\delta_{x}$ denotes the Dirac delta function centered at $x$. Thus

[TABLE]

which corresponds to the kinetic energy of the particle integrated over time. Likewise one computes

[TABLE]

where $d_{i,j}$ denotes the minimal distance between the point $X(t_{i})$ and the LOR $l_{j}$.

Let us for a moment ignore the contribution of scattering, i.e. we set $p=0$. If the pair $(\rho,\omega)$ corresponding to a single particle is inserted into (5) one obtains

[TABLE]

For fixed $m$ the data term in this functional is (up to a factor) the sum of the squared distances $d_{i,j}^{2}$ just as in [6]. The regularization term $\int_{0}^{T}\|X^{\prime}(t)\|_{2}^{2}\,{\mathrm{d}}t$ leads to the preference of short trajectories with low derivative. While the precise term used in [6] is different (it penalizes quadratically the change in the cubic B-spline coefficients in a parametrization of $X$) it is very similar in spirit and can in fact be interpreted as a time discretization of $\int_{0}^{T}\|X^{\prime}(t)\|_{2}^{2}\,{\mathrm{d}}t$.

Now consider scattering. One finds $A^{\mathrm{s}}_{i,j}\rho={\mathrm{const}}\cdot m$, thus

[TABLE]

with $f(d)=-\log(\exp(\tfrac{-d^{2}}{2\epsilon^{2}})+\tilde{p})$ where $\tilde{p}\geq 0$ is the product of $p$ and some constant. For $d/\epsilon\ll-\log\tilde{p}$ one has $f(d)\approx\tfrac{d^{2}}{2\epsilon^{2}}$, while for $(d/\epsilon)^{2}\gg-\log\tilde{p}$ one finds $f(d)\approx-\log\tilde{p}$. Thus, the function $g(d)=\min\{\tfrac{d^{2}}{2\epsilon^{2}},-\log\tilde{p}\}$ can be viewed as a simple approximation of $f$, and indeed the function $g$ is used in [6] to model scattering and random coincidences.

Summarizing, when restricted to a single particle our functional is indeed very similar to the one proposed in [6]. Analogously to the above computations one could now derive a generalization of their functional for multiple particles with trajectories given by functions $X_{k}$ and masses $m_{k}$ for $k=1,\ldots,P$. However, the resulting functional would be severely non-convex in the coordinates $X_{k}$ and thus difficult to optimize. The main novelty of our functional is thus not the forward model or the regularization term, which are similar to the ones used in [6], but indeed the formulation in terms of time-dependent densities and flux fields, which yields a convex functional that is readily minimized and does not require to manually choose the particle number $P$.

III Methods

III-A Discretization and numerical optimization

The functions $\rho_{t}$ and $\omega_{t}$ are discretized by a regular grid of voxels over $\Omega$ for time points $t\in\{0,\Delta T,\ldots,T\}$. Constraint (2) and $\mathcal{S}(\rho,\omega)$ are then discretized using finite differences on a staggered grid as in [13]. The operators $A_{i,j}^{\mathrm{d}}$ and $A_{i,j}^{\mathrm{s}}$ are discretized by identifying for each pair $j$ of detectors which voxels contribute how much. In particular the discretization of the reconstruction method is completely independent from the software used to generate the simulated data, thus preventing the fundamental inverse crime [14]. Therefore, the obtained results should indeed be meaningful.

Unless stated otherwise, in our experiments we set $T=120\,\textnormal{s}$, $\Omega=[0,L]\times[0,L]\times[0,L/4]$ with $L=160\,\textnormal{mm}$ and discretize $[0,T]\times\Omega$ with a regular Cartesian grid with $65\times 64\times 64\times 16$ points, resulting in a meshsize of $\Delta L=2.5\,\textnormal{mm}$. We set the width $\epsilon$ of the Gaussian kernel $G(y,x)$ from (1) to $\epsilon=5\textnormal{mm}$ which is slightly larger than the value in the simulations (cf. fig. 2) to reduce discretization artifacts.

The functional $J^{E,p}$ to be minimized is jointly convex in $(\rho,\omega)$ (indeed, the constraints $\rho\geq 0$ and (2) are linear, and $\mathcal{S}$ is known to be convex as well [12]). Thus, its minimization can be performed via a convex optimization approach. The optimization problem is then written in the form

[TABLE]

where $K$ is a suitable matrix and $F$ and $H$ are convex lower-semicontinuous functionals. To this optimization problem the primal-dual implicit gradient descent and ascent from [15] is applied. The precise form of $F,H$ and $K$ as well as the algorithm are given in the supplementary materials.

III-B Framewise reconstruction

To demonstrate the benefit of temporal regularization we implement framewise reconstruction as a reference, i. e. we set $\beta\!=\!0$ and assume $\rho$ to be constant on each interval $\tau_{i}$. The flow field $\omega$ no longer appears in the functional. It is trivial to adapt the numerical scheme accordingly. Apart from our scatter-unbiasing the unregularized functional is then comparable to standard Bayesian reconstruction methods for PET imaging.

We expect that for small $\Delta T$ the framewise method suffers from a lack of signal in each independent time frame. As $\Delta T$ increases, more events become available in each frame, but there will be increasing inaccuracy due to particle motion.

III-C Quantification of reconstruction error

Using simulated measurements allows us to quantitatively compare our reconstruction results to the ground truth. We have to simultaneously measure errors in the localization of multiple particles and their respective masses (since PET is quantitative). Note that both cannot be separated from each other, since incorrectly positioning some mass from $x\in\Omega$ at $y\in\Omega$ can be interpreted both as an error in localizing the mass from $x$ or as an error in reconstructing the correct amount of mass in $x$ and $y$. In particular we cannot use the squared Euclidean distance from [6] as it only applies to single particles and cannot take into account masses. A suitable joint error quantification for particle localization and mass is achieved by the Wasserstein–Fisher–Rao (WFR) metric $d_{{\mathrm{WFR}},\alpha}(\rho,\tilde{\rho})$ between material distributions $\rho$ and $\tilde{\rho}$ (see supplementary material or [11] for a precise definition). It measures how far on average the material in $\rho$ has to be transported to change $\rho$ into material distribution $\tilde{\rho}$, where transport beyond $\text{distance}\approx\alpha$ is replaced by just changing the mass locally (assuming that rather the mass intensity was wrong than its localization). In detail, for two Dirac masses $\delta_{x}$ at $x$ and $\delta_{y}$ at $y$ one finds

[TABLE]

where $\alpha$ is a length scale parameter that balances the trade-off between localization and mass errors. This can be interpreted as a truncated Euclidean distance since $d_{{\mathrm{WFR}},\alpha}^{2}(\delta_{x},\delta_{y})\approx\|x-y\|^{2}$ for $\|x-y\|\ll\alpha$. Consequently, it provides a natural generalization of the Euclidean distance to the setting of multiple particles with different masses.

To quantify our reconstruction errors we pick $\alpha=25\,\textnormal{mm}$ and calculate

[TABLE]

for $\rho^{\mathrm{gt}}$ the ground truth, $\rho$ our reconstruction, and $\mathcal{N}$ the map that normalizes a non-negative density. We apply normalization to avoid tedious calibration of the absolute scanner sensitivity. ${\mathrm{err}}$ can roughly be interpreted as mean localization error for the particles with a truncation for deviations $\gg\alpha$. It is computed numerically using the code from [16].

Since $\alpha\gg\Delta L$ (where $\Delta L$ is the distance between two neighbouring grid points) the error inflicted by misplacing all particles by one discrete grid point is approximately

[TABLE]

This provides a scale for the expected discretization error.

IV Experimental validation

IV-A Monte Carlo Simulation

We work on simulated data. In this way we can tune the problem difficulty, and reliable ground truth is available, allowing a transparent evaluation of our reconstruction method.

We simulate list-mode data of a clinical PET scanner with the Geant4 Application for Tomographic Emission (GATE) Monte Carlo tool [17]. Within this simulation, the geometry and detector concept of the Siemens Biograph mCT PET/CT is implemented to mimic a state-of-the-art human PET system with LSO detector material [18]. To simulate cells under realistic conditions, body anatomy and attenuation is taken from a 4D extended cardiac-torso (XCAT) phantom [19]. A single cell is simulated as a point source (diameter 1 mm) with radioactivity of 10 kBq 18F (${T_{\frac{1}{2}}}=6586\,\textnormal{s}$). This cell is moved in 120 steps along a circular trajectory with radius 60 mm in the x-y plane, see fig. 1. At each step, a Monte-Carlo simulation is performed using 1 s simulation time, and the information of the detected coincidence events is stored in a list-mode data format. Note that due to this discretization of a continuous motion the simulated cell moves abruptly in steps of $\pi\,\mathrm{mm}\approx 3.14$ mm which will cause a slight additional error in our reconstructions (so we actually overestimate our reconstruction error in the experiments).

Two different scenarios are simulated. In the first scenario (‘without scattering’), $\beta^{+}$ particles are simulated with typical energies of 18F followed by positron-electron annihilation inside the human tissue and emission of two 511 keV gammas. Scattered gamma events are rejected by setting an energy window of 510–650 keV with no energy blurring. For a 10 kBq source activity the rate of counts detected by the scanner was $44\,\textnormal{cps}$ (counts per second), corresponding to a ratio of $s=4.4\cdot 10^{-3}\,\textnormal{cps}/\textnormal{Bq}$. In the second scenario (‘with scattering’), scattered events are partly allowed using an energy window of 435–650 keV (energy blurring 11.6 % at reference energy 511 keV), a range which is typically used in a clinical setup. The intrinsic natural radiation from Lu-176 within the LSO crystal was not considered in the GATE simulation. Events were detected with a rate of $104\,\textnormal{cps}$, yielding a ratio of $s=1.0\cdot 10^{-2}\,\textnormal{cps}/\textnormal{Bq}$. The scatter fraction in this scenario was about 18%. In both simulations a coincidence window of 4.1 ns is used.

To obtain more realistic experimentally feasible conditions the activity is then subsampled to rates between 160 Bq and 800 Bq in the first simulation scenario (corresponding to 0.7 to 3.6 cps for a single particle) and 500 Bq in the second scenario (5.2 cps/particle). Multiple particles are simulated by applying a rotation in the x-y plane to the listmode data and then combining the resulting events. This is admissible since the scanner geometry is approximately symmetric around this axis. Different subsamplings of the original data were used so that events from different particles are statistically independent.

To obtain a notion of robustness of our reconstruction method we generated five runs for each simulation scenario. The reconstruction accuracies and scattering ratios reported in the following are the obtained averages, and we visualize the standard deviation among the five runs with errorbars.

IV-B Basic setup

Our primary test phantom consists of four cells moving one after another on a circle in the x-y plane with a radius of 60 mm. We consider four experimental parameters: the amount $m$ of radionuclides within the cells, the cell velocity $v$, the distance between subsequent cells and whether or not photon scattering is included in the recorded events. The detection rate (counts per second) is then given by $A=s\cdot m\cdot\log 2/{T_{\frac{1}{2}}}$. We vary $A$ by varying $m$. For values of $s$ and ${T_{\frac{1}{2}}}$ and how $m$ is varied by subsampling the original data, see section IV-A.

Recall that in addition there are two major parameters in the reconstruction model: $\beta$ for kinetic regularization and the scattering tuning parameter $p$. This setup allows to carefully test the method at different levels of difficulty and to understand how to choose the two reconstruction parameters properly.

Figure 2 shows the cumulative distributions of the distance from the original positron emission to the line of response of the activated detector pair for all detected events, for simulated data with and without photon scattering. In absence of scatter this distance is caused by the positron range and photon non-collinearity, as modelled with the kernel $G$ in $A^{\mathrm{d}}_{i,j}$, see (1). With scatter included, minor scattering with small perturbations to the photon orientation is still modelled with $G$. In events with strong scattering the true photon path will differ substantially from the idealized line of response, and thus the distance between positron emission and line may be much larger. This is modelled with the operator $A^{\mathrm{s}}_{i,j}$, and the long tail of the scattered distribution indicates that the assumption of ‘uniform’ scattering is indeed a reasonable approximation.

IV-C Temporal regularization

We now study the influence of detection rate and temporal regularization. For simplicity we start without scattering.

Exemplary visualizations of reconstructed trajectories with and without temporal regularization are given in figs. 3 and 4. As the particles in the phantom are restricted to the x-y plane, for simplicity we visualize the projection of $\rho$ into this plane (or equivalently, the integration of $\rho$ along the z axis).

Figure 5 shows the reconstruction error for different detection rates $A$ and values of $\beta$, where the detection rate is simply varied by changing the amount of radioactive material in the simulation. The optimal value of $\beta$ is found to lie between 0.4 s/mm2 and 1.0 s/mm2 and is essentially independent of the rate $A$ and number of time frames $M$. For detection rates of $A=5.7\,\textnormal{cps}$ and above the optimal reconstruction error is almost constant and almost at the level of discretization artifacts (see (6)) across about one order of magnitude for $\beta$. For lower rates the error increases more quickly. The optimal reconstruction error is approximately proportional to the inverse activity (cf. fig. 7 (right)). For $A=(2.8,4.3,8.3)\,\textnormal{cps}$ the lowest reconstruction errors are approximately $(8.0\pm 0.7,5.3\pm 0.5,3.6\pm 0.1)\,\textnormal{mm}$. For low regularization $\beta$ the error increases since the coupling between the information contained in different time frames becomes weaker (this effect is more pronounced with more time frames, $M=129$, each of which then contains less information). As $\beta$ increases the error grows due to overregularization. Fluctuation between different simulation runs is reasonably small (as the errorbars indicate) showing that our method is robust.

For comparison, the reconstruction errors without kinetic regularization (i.e. framewise reconstruction, section III-B) are given in fig. 6. The errors are higher than with temporal regularization, and in particular they increase drastically with decreasing detection rate and therefore also with increasing number of time frames (which reduces the available information per frame). Reducing the number of frames increases the available information per frame but the error eventually increases due to ‘motion blur’. The reconstruction with temporal regularization is more robust to little information per frame due to the coupling between the frames.

Parameter $\beta$ balances the strength of temporal regularization versus agreement with the measurements. Its optimal value depends on the substance amount $m$, particle velocity $v$ and also the halflife ${T_{\frac{1}{2}}}$ (since our data term uses radioactivity, but the regularization is based on material mass and thus still needs to be scaled by the inverse halflife). The precise relation can be deduced from the scaling behavior of the functional $J^{E,p}$.

For instance, if one rescales $m\to q\cdot m$, ${T_{\frac{1}{2}}}\to q\cdot{T_{\frac{1}{2}}}$ for some $q>0$ (which leaves the rate $A$ and thus the measurement $E$ unchanged) and rescales $\beta\to\beta/q$, then the minimizer of $J^{E,p}$ is rescaled according to $(\rho,\omega)\to(q\cdot\rho,q\cdot\omega)$, i. e. the optimal reconstruction is rescaled by the same factor as the true material amount. As a consequence the reconstruction quality is invariant under this rescaling and thus if $\beta$ is the optimal parameter for $(m,v,{T_{\frac{1}{2}}})$ then $\beta/q$ is the optimal parameter for $(q\cdot m,v,q\cdot{T_{\frac{1}{2}}})$. With additional similar arguments one can show that the optimal parameter $\beta$ is given by $C/({T_{\frac{1}{2}}}\cdot v^{2})$ for a suitable constant $C$. From fig. 5 we deduce that the optimal $\beta$ for $v=3.14\,\textnormal{mm}/\textnormal{s}$ and ${T_{\frac{1}{2}}}=6586\,\textnormal{s}$ is approximately $0.9\,\textnormal{s}/\textnormal{mm}^{2}$ and therefore $C\approx 600$. Unless stated otherwise, in all subsequent experiments we set $\beta$ according to this formula.

In summary, one obtains good reconstructions for a large interval of values of $\beta$, and the dependence of the optimal value on the experimental parameters is well understood. Also, as compared to the framewise reconstruction, the regularized reconstruction quality is more robust at low detection rates.

IV-D Particle velocity and tracking efficiency

In [6] it was observed that particle velocity $v$ and detection rate $A$ are not independent parameters with respect to the reconstruction quality. In essence only their ratio $v/A$ is important which was called tracking efficiency in [6]. This is intuitive as the tracking efficiency is the average distance between two detected positron decays along the particle trajectory.

This relation can also be observed for our reconstruction functional $J^{E,p}$. Let $q>0$ and rescale $v\to q\cdot v$ and $m\to q\cdot m$ which implies $A\to A\cdot q$ and in particular that $v/A$ remains constant. Then, similar to above, if one rescales $\beta\to\beta/q^{2}$, one can show that the minimizer of the functional transforms as $\rho_{t}\to q\cdot\rho_{q\,t}$ and in particular the reconstruction error does not change. This is confirmed numerically in fig. 7 (left). As predicted the curves look very similar.

The relation between tracking efficiency and tracking accuracy for the data of fig. 5 is shown in fig. 7 (right). For a detection rate of $A=4.3\,\textnormal{cps}$ for four particles the detection rate per particle is $A_{\textnormal{single}}=1.1\,\textnormal{cps}$ which corresponds to a tracking efficiency of $3.0\,\textnormal{mm}$, at a tracking accuracy of $5.3\,\textnormal{mm}$. In [6] comparable tracking efficiencies are reported at somewhat better accuracies albeit only for a single particle and a smaller scanner geometry providing higher image resolution.

IV-E Scattering

We now turn to data with scattering. We will use the full forward operator involving the parameter $p$ used for unbiasing.

Figure 8 (left) shows the reconstruction error for two values of $p$. For low and medium values of $\beta$ the behavior is very similar to the no-scattering case (cf. fig. 5). However, with scattering the effect of overregularization is more dramatic which is due to a new mechanism that becomes available with scatter modelling: when the kinetic regularization is too strong, the lowest functional values can be obtained by keeping the reconstructed particles almost stationary and declaring almost all detected events as scattering. Figure 8 (right) shows the corresponding estimated scattering ratios, where we interpret an event $(i,j)$ as ‘scattered’ if $pA_{i,j}^{\mathrm{s}}\rho\geq A_{i,j}^{\mathrm{d}}\rho$. For small and medium values of $\beta$ the ratio stays almost constant for each $p$, but increases quickly to 1 in the overregularization regime.

For sufficiently low $\beta$ (below overregularization) the scattering ratios for $p=0.001$ and $p=0.002$ are roughly 0.2 and 0.8, respectively. From fig. 2 we see that the former is relatively accurate while the latter is substantially too high and wrongfully discards many unscattered events. Despite this, the reconstruction error for $p=0.002$ is almost on par with $p=0.001$ since the remaining 20 % of undiscarded events provide sufficient information about the particle locations.

Due to the stronger effect of overregularization the optimal value for $\beta$ is slightly reduced. Hence, for scattered data we will set in the following $\beta=C/({T_{\frac{1}{2}}}\cdot v^{2})$ for $C\approx 130$.

Figure 9 illustrates the influence of $p$ on the reconstruction error and estimated scattering ratio in more detail. The estimated scattering ratio increases monotonically with $p$. For small $p$ only the ‘most absurd’ events are labeled as scattered. With $p$ the ratio increases, and eventually almost all events are considered scattered. In accordance, for small $p$ the reconstruction error is higher since some scattered events that are not correctly recognized as such cause artifacts in the reconstruction. For high $p$ the error increases since too many events are discarded as scattered which leaves less reliable information for reconstruction. The minimum reconstruction error is roughly obtained for that value of $p$ which predicts the correct scatter ratio. Our experiments show that the optimal value for $p$ does not depend significantly on detection rate or particle velocity, and it is thus relatively easy to calibrate. If $p^{\mathrm{d}}$ and $p^{\mathrm{s}}$ are known, the nuclide intensity can in principle be reconstructed accurately, but regularization adds a systematic bias towards lower masses. For the chosen optimal values of $\beta$, the reconstructed mass was around 87% of the true value.

Reconstructions of data with scattering, with and without proper detection of scattered events are visualized in fig. 10.

IV-F Particle distance

In fig. 11 the reconstruction errors for different detection rates with and without scattering are investigated for varying distances between the particles. For sufficiently large distances the reconstruction error is low and approximately constant, indicating that particles can be separated cleanly. As particles come closer the reconstruction error increases. Eventually it decreases again when the particles are essentially merged.

Visualizations of reconstructions for various particle distances for the data of fig. 11 with scattering are shown in fig. 12.

V Conclusion

We presented a mathematical model and algorithm for the reconstruction of a temporally changing radionuclide distribution within a patient from PET listmode data. As opposed to classical clinical PET reconstruction, our method performs a dynamic rather than static imaging, resolving the temporal changes. It is in particular designed for tracking individual radiolabelled cells, but could without modification also be applied to the reconstruction of a temporally evolving distribution of injected PET tracers or PET reconstruction under organ (cardiac or respiratory) motion.

In essence, for tracking individual cells via PET there are two existing competing approaches, time binning with independent reconstruction for each bin and path fitting [6]. Time binning pretends the radionuclide distribution was stationary over short time intervals and performs a classical reconstruction over those time intervals. As a drawback, too large time intervals cause substantial blurring due to the radionuclide motion, while too small time intervals contain too few detected decays to allow a reliable image reconstruction separately for each time interval. In contrast, our novel method is essentially continuous in time so that no blurring occurs, but at the same time exploits temporal consistency so that the reconstruction almost behaves as if the radionuclide distribution was indeed stationary. The resulting strong superiority over time binning was demonstrated quantitatively on simulated data.

Path fitting on the other hand tries to fit a particle trajectory to the detected decay events. While for a single particle this is readily doable, for multiple particles it introduces the combinatorial problem of deciding which event belongs to which particle. Our approach can be seen as a relaxation of this method in the sense that in principle our method also allows particles to split or to merge. As a consequence our reconstruction boils down to a simple convex optimization.

The above advantages come at the price of two model parameters that have to be adequately picked, however, the influence of either parameter is well-understood. One parameter determines the amount of detected events which are interpreted as having undergone scattering (or coming from background noise) and are thus discarded. Since the amount of scattering and noise is known for each setting (it depends on the PET scanner, the radionuclide and the imaged object), this parameter can be easily tuned (once and for all for each setting) by comparing the estimated scatter ratio to the expected one. The other parameter represents the regularization strength, and we have shown how to choose it optimally depending on radionuclide halflife and expected particle velocity.

We tested our method on simulated data for imaging a human in a standard clinical PET scanner, using four cells or particles, each with 160-800 Bq activity and 3.1 mm/s velocity (reconstructions of larger numbers of cells would exhibit the same quantitative behaviour). For example, at 240 Bq per cell (corresponding to 1.1 cps per cell), four cells could be tracked simultaneously with an accuracy of 5.3 mm. The particle trajectories can clearly be identified and separated up to a distance comparable to the typical positron range (roughly two to five times, depending on the activity).

The chosen cell test setting is highly challenging but not completely unrealistic under the assumption that a single cell can be labelled with an activity of $\sim 200\,\textnormal{Bq}$ per cell. Although the activity is usually in the range of a few Bq per cell for standard labelling protocols (e.g. [20, 21]), it was already demonstrated that high activity levels in the range of 50-70 Bq/cell could be achieved using F-18-FDG or Ga-68 labelled cancer cells [10, 22]. PET isotopes with short half-life would be preferable for tracking fast moving cells as less molecules need to be attached to the cell providing high activity. Conversely, long living isotopes, such as Zr-89 with a half-life of 78.4 h, would allow to track cells over several days albeit with limited activity per cell. Thus, the choice of isotope and labelling efficiency will be crucial for tracking accuracy. Furthermore, moving from a simulated clinical PET scanner to a total-body human scanner, of which prototypes are under evaluation, the sensitivity and thus detection rate could be increased substantially (a 40-fold effective sensitivity was predicted in [23]). This means that the radioactivity per cell could be reduced accordingly or particles could be tracked at higher speed, e.g. also inside large vessels such as the aorta.

Finally, reconstruction results could be further improved by employing more a priori knowledge. Here we completely ignored any anatomic information (which could come from PET-MR or PET-CT systems), and rather than regularizing the kinetic energy quantities such as the acceleration might be biologically more relevant. Also, any time-of-flight (TOF) information of the coincidence events is ignored in our model and its incorporation would require weaker regularization as the reconstruction problem becomes less ill-posed and the signal-to-noise characteristics are substantially improved (as is typically observed in modern PET systems). Thus, the bias introduced by our transport regularization will be reduced, enhancing the tracking accuracy of our algorithm further.

Supplementary material

Numerical optimization

We will perform numerical optimization of the discretized problem with the primal-dual implicit gradient descent and ascent from [15]. The (primal) variables of the problem are $\rho$ and $\omega$, the mass and flows on the time- and space-staggered grids (cf. [13]) which we summarize as $x=(\rho,\omega)^{\top}$. The discrete problem can then be written as

[TABLE]

where $K$ is a combination of linear operators given by

[TABLE]

Note that $F$ and $H$ are convex and lower-semicontinuous. The dual variables corresponding to the arguments of $F$ will be denoted by $y=(\varphi_{1,\rho},\varphi_{1,\omega},\varphi_{2},\varphi_{3})$. Let $\tau,\sigma>0$ with $\tau\cdot\sigma\|K\|^{2}<1$ where $\|K\|$ denotes the operator norm of $K$. Then, for initial primal and dual variables $x^{(0)}$ and $y^{(0)}$ the algorithm is given for $\ell=0,1,\ldots$ by

[TABLE]

with the convention $x^{(-1)}=x^{(0)}$. Here $F^{\ast}$ denotes the Fenchel–Legendre conjugate of $F$ and $\mathrm{Prox}_{\sigma F^{\ast}}$ denotes the proximal operator of $F^{\ast}$ with stepsize $\sigma$, defined by

[TABLE]

and analogously for $H$. The Fenchel–Legendre conjugate of $F$ decomposes into the conjugates of $F_{1}$, $F_{2}$ and $F_{3}$,

[TABLE]

and so does the proximal operator,

[TABLE]

Each of these proximal operators is separable into independent operators per grid point or detector bin. The proximal operator $F^{\ast}_{1}$ is described in [13], those for $F^{\ast}_{2}$ and $F^{\ast}_{3}$ can be determined by elementary calculations. The proximal operator of $H$ amounts to computing the projection onto solutions of the continuity equation which involves the solution of the Poisson equation on a Cartesian grid. This can be performed efficiently with the FFT algorithm [13].

Wasserstein–Fisher–Rao metric

The WFR metric between two nonnegative densities $\mu,\nu:\Omega\to[0,\infty)$ is defined via

[TABLE]

Above, the optimization variables are the temporally changing material density $\rho_{t}$, the temporally changing vector field $\omega_{t}$ describing a material transport, and $\zeta_{t}$ describing a temporal change of material mass. For simplicity, in our notation we treated $\rho_{t}$, $\omega_{t}$, and $\zeta_{t}$ as functions, but in the mathematically rigorous definition the minimization is performed in the space of Radon measures, i.e. $\rho$, $\omega$, and $\zeta$ are all measures on $[0,T]\times\Omega$ (for details see [11]).

By definition $d_{{\mathrm{WFR}},\alpha}^{2}(\mu,\nu)$ is the minimum action to move $\mu$ onto $\nu$, where in contrast to $\mathcal{S}$ the action now is complemented with a term for mass changes during transport. The parameter $\alpha>0$ controls how much mass change is penalized compared to transport. It can be interpreted as the length scale across which a mass misplacement can be interpreted as a localization error (in fact, it is known that mass transport only happens up to a distance $\pi\cdot\alpha$ in $d_{{\mathrm{WFR}},\alpha}$).

For two Dirac measures at $x$ and $y$ with masses $m_{x}$ and $m_{y}$ one obtains

[TABLE]