A probabilistic framework for uncertainty quantification in positron emission particle tracking

The space of lines in $\mathbb{R}^3$ is 4-dimensional since lines can be characterized by a unit vector n , giving their direction, and the closest point to the origin, which lies in the plane normal to n . With the restriction $n_z \geqslant 0$, this identifies the space of lines with the tangent bundle to a hemisphere. We parameterize n using the standard polar spherical angles. The point on the line closest to the origin is then conveniently represented as

$$\begin{equation} \boldsymbol{a}=a_{\varphi}\boldsymbol{e}_{\varphi}+a_{\theta}\boldsymbol{e}_{\theta} \end{equation} \tag{ 1 }$$

in terms of the unit vectors $\boldsymbol{e}_{\varphi}$ and $\boldsymbol{e}_{\theta}$ of the spherical coordinate system (with $\boldsymbol{e}_r = \boldsymbol{n}$, see figure 2). In this way, a LoR L is parameterized as

$$\begin{equation} L=\left(\varphi,\theta,a_{\varphi},a_{\theta}\right)\in\left[0,2\pi\right]\times\left[0,\pi/2\right]\times\mathbb{R}^{2}. \end{equation} \tag{ 2 }$$

There is a natural measure on that space, $\mathrm{d} \mu$ say, characterized by invariance in rotation and translation [16]. In terms of the coordinates (2), this has the simple form

$$\begin{equation} \mathrm{d} \mu = \sin \theta \, \mathrm{d} \varphi \mathrm{d} \theta \mathrm{d} a_\varphi \mathrm{d} a_\theta. \end{equation} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

The generation of LoRs given the positions and activities of K radioactive particles defines a Poisson point process in the space of lines and in time. This has long been known in PET, with numerous works describing LoR generation as a Poisson process (see, e.g. [17, 18]), however in PEPT the radioactive source position is inherently time dependent. We denote by $\lambda(L|\boldsymbol{x}(t),\rho)$ the rate (or intensity) of the Poisson process corresponding to the generation of LoRs by a single particle with activity ρ moving along the path $\boldsymbol{x}(t)$. Thus $\lambda \mathrm{d} \mu \mathrm{d} t$ is the probability of finding a LoR in the volume $\left[\varphi,\varphi+\mathrm{d}\varphi\right]\times\left[\theta,\theta+\mathrm{d}\theta\right]\times\left[a_{\varphi},a_{\varphi}+\mathrm{d} a_{\varphi}\right]\times\left[a_{\theta},a_{\theta}+\mathrm{d} a_{\theta}\right]$ during a time interval $\left[t,t+\mathrm{d} t\right]$. The form of $\lambda(L|\boldsymbol{x}(\cdot),\rho)$ depends on the geometry of the detectors and physics of positron annihilation, scattering and other processes; we give an explicit form for a cylindrical detector in section 3. Note that $\lambda(L|\boldsymbol{x}(t),\rho)$ depends on time implicitly through the time dependence of the particle position x . The rate of generation by K particles is simply $\sum_{k = 1}^K$ $\lambda(L|\boldsymbol{x}_k(t),\rho_k)$.

Using the standard formula for Poisson point processes [19], we can write down the probability of observing N LoRs $L_1,\ldots,L_N$ at times $0\lt t_1\lt\cdots\lt t_N\lt\Delta t$ as

$$\begin{equation} P\left(\mathcal{L}_t|\boldsymbol{X}(\cdot),\varrho,\Delta t\right)=\mathrm{e}^{-\Lambda(\boldsymbol{X}(\cdot),\varrho,\Delta t)}\prod_{n=1}^N \sum_{k=1}^K \lambda\left(L_{n} |\boldsymbol{x}_k(t_n),\rho_k \right). \end{equation} \tag{ 4 }$$

Here $\mathcal{L}_t = \left\{(L_{1},t_1),\ldots,(L_{N},t_N)\right\}$ denotes the observed LoRs and times of observation and we have defined

$$\begin{equation} \Lambda(\boldsymbol{X}(\cdot),\varrho,\Delta t)= \sum_{k=1}^K \int_0^{\Delta t} \int \lambda\left(\left.L\right|\boldsymbol{x}_k(t),\rho_k\right)\mathrm{d}\mu \mathrm{d} t, \end{equation} \tag{ 5 }$$

with the integration with respect to $\mathrm{d} \mu$ carried out over the entire space of LoRs. Note that $P\left(\mathcal{L}_t|\boldsymbol{X}(\cdot),\varrho,\Delta t \right)$ is properly normalized, in the sense that

$$\begin{equation} \sum_{N=0}^\infty \int_0^{\Delta t} \mathrm{d} t_1 \int_{t_1}^{\Delta t} \mathrm{d} t_2 \cdots \int_{t_{N-1}}^{\Delta t} \mathrm{d} t_N \int\cdots\int P\left(\mathcal{L}_t|\boldsymbol{X}(\cdot),\varrho,\Delta t \right) \, \mathrm{d} \mu_1 \cdots \mathrm{d} \mu_N =1. \end{equation} \tag{ 6 }$$

In most inversion algorithms, the time interval $\Delta t$ is taken small enough that particles can be assumed fixed (see [20], however). With this assumption, information about the times of observations of the LoRs can be disregarded, and we can focus on the probability of observing a set $\mathcal{L} = \{L_1,\ldots,L_N\}$ of LoRs, irrespective of the times. This probability is obtained by integrating (4) with respect to $0\lt t_1 \lt t_2 \lt \cdots \lt t_N \lt \Delta t$ to find

$$\begin{equation} P\left(\mathcal{L}|\boldsymbol{X},\varrho,\Delta t\right)=\frac{(\Delta t)^N \mathrm{e}^{-\Lambda(\boldsymbol{X},\varrho,\Delta t)}}{N!} \prod_{n=1}^N \sum_{k=1}^K \lambda\left(L_{n} |\boldsymbol{x}_k,\rho_k \right), \end{equation} \tag{ 7 }$$

where the notation emphasizes that X and x _k are regarded as constant. Equation (5) also simplifies to

$$\begin{equation} \Lambda(\boldsymbol{X},\varrho,\Delta t)=\Delta t \sum_{k=1}^K \int \lambda\left(\left.L\right|\boldsymbol{x}_k,\rho_k\right)\mathrm{d}\mu. \end{equation} \tag{ 8 }$$

Bayes' formula provides the means to infer a probability distribution for the particle positions and activities from LoR data. Under the assumption of fixed particles, it takes the form

$$\begin{equation} P(\boldsymbol{X},\varrho |\mathcal{L},\Delta t) \propto P(\mathcal{L}|\boldsymbol{X},\varrho,\Delta t)P(\boldsymbol{X},\varrho), \end{equation} \tag{ 9 }$$

where $P(\boldsymbol{X},\varrho |\mathcal{L},\Delta t)$ is the posterior probability for the (fixed) particle positions $\boldsymbol{X} = \left\{\boldsymbol{x}_1,\ldots,\boldsymbol{x}_K\right\}$ and activities $\varrho = \left\{\rho_1,\ldots,\rho_K\right\}$, $P(\boldsymbol{X},\varrho |\mathcal{L},\Delta t)$ is the likelihood given in (7), and $P(\boldsymbol{X},\varrho)$ is the prior. In the application of section 4, we take a uniform (improper) prior over the set of physically realizable particle positions and activities, namely

$$\begin{equation} P\left(\boldsymbol{X},\varrho\right)= \begin{cases} 1 & \textrm{if} \ x_{k}\in V \ \textrm{and} \ \rho_{k}\geqslant 0 \quad \forall k \\ 0& \textrm{otherwise} \end{cases}, \end{equation} \tag{ 10 }$$

where V is the volume enclosed by the detectors.

The assumption of fixed particle positions that underlies (9) restricts the length of the observation time $\Delta t$ and hence of the number of observed LoRs. Depending on the speed of the particles, this can severely limit the accuracy of the inversion. To overcome this, we can retain information about the times of LoR observations and carry out Bayesian inference for the particle trajectories rather than fixed positions, using (4) as the likelihood. In practice, the trajectories must be parameterized, say as $\boldsymbol{x}_k(t) = \hat{\boldsymbol{x}}(t,\vartheta_k)$, where ϑ_k groups the parameters of the trajectory of particle k. Bayes' formula is then used in the form

$$\begin{equation} P(\boldsymbol{\vartheta},\varrho |\mathcal{L}_t,\Delta t) \propto P(\mathcal{L}_t|\boldsymbol{\vartheta},\varrho,\Delta t)P(\boldsymbol{\vartheta},\varrho), \end{equation} \tag{ 11 }$$

where $\boldsymbol{\vartheta} = \{\vartheta_1,\ldots,\vartheta_k\}$ and the likelihood $P(\mathcal{L}_t|\boldsymbol{\vartheta},\varrho,\Delta t)$ is obtained from (4) by substituting the parametric form of the trajectories. A simple parameterization can be constructed by a Taylor expansion of the particle trajectory around the middle of the observation time. We follow Manger et al [9] by writing

$$\begin{equation} \boldsymbol{x}_k\left(t\right)=\sum_{m=0}^{M} \frac{1}{m!} \left.\frac{\mathrm{d}^{m}\boldsymbol{x}_k}{\mathrm{d}t^{m}}\right|_{\Delta t/2}\left(t-\Delta t/2\right)^{m}+O\left(\Delta t^{M+1}\right). \end{equation} \tag{ 12 }$$

and use the position and its M first derivatives as parameters:

$$\begin{equation} \vartheta_k=\left( \left. \boldsymbol{x}_k \right|_{\Delta t/2},\ldots,\left. \frac{\mathrm{d}^M\boldsymbol{x}_k}{\mathrm{d} t^M} \right|_{\Delta t/2} \right). \end{equation} \tag{ 13 }$$

We refer to inference of this type as high-order inference and to M as the order of the inference. We implement the M = 1, first-order inference in section 4.4.

The LoRs are detected on a surface of discrete detectors, in our case a cylindrical envelope. The line through x intersects this cylinder at two points which we parameterize as

$$\begin{equation} u^{^{\prime}}=\left(\phi^{^{\prime}}_{1},\phi^{^{\prime}}_{2},z^{^{\prime}}_{1},z^{^{\prime}}_{2}\right)\in\left[0,2\pi\right]\times\left[0,2\pi\right]\times\left[-H/2,H/2\right]\times\left[-H/2,H/2\right], \end{equation} \tag{ 15 }$$

where $\phi^{^{\prime}}_{1,2}$ and $z^{^{\prime}}_{1,2}$ are the azimuthal angles and heights of the two points (see figure 1 for illustration of this parameterization). Measurement uncertainty introduces a deviation between the detected coordinates u of these intersections and the exact intersection uʹ. As a result, the coordinates $L = (\varphi,\theta,a_\varphi,a_\theta)$ of the observed LoR differ from the coordinates $L^{^{\prime}} = (\varphi^{^{\prime}},\theta^{^{\prime}},x_{\varphi^{^{\prime}}},x_{\theta^{^{\prime}}})$ of the line followed by the photons. The uncertainty associated with the finite size of the detectors implies that the intersection parameters u of the observed LoR are distributed uniformly over the area of the two detectors intersected by the true LoR. To avoid discontinuous distribution and to group all measurement errors in a single, simple model, we replace this uniform distribution by a normal distribution: we assume that u is normally distributed about uʹ according to

$$\begin{equation} u | L^{^{\prime}} \sim \mathcal{N}(u^{^{\prime}},C), \quad \textrm{where} \ \ C=\left(\begin{array}{cccc} g^{2}/R^{2} & 0 & 0 & 0\\ 0 & g^{2}/R^{2} & 0 & 0\\ 0 & 0 & g^{2} & 0\\ 0 & 0 & 0 & g^{2} \end{array}\right). \end{equation} \tag{ 16 }$$

Here g² is the variance of the detected heights z₁ and z₂; assuming that the uncertainty about the true LoR is isotropic, it follows that the variance in $\phi_{1},\phi_{2}$ is $g^2/R^2$. In practice, we fix g to ensure the Gaussian distribution does not strongly deviate from a top hat function; specifically, we take g such that the difference between integrating the top hat and Gaussian distributions over a single detector is minimized.

There is a one-to-one map between the coordinates of LoRs and the coordinates u of the detected points which we denote by $u = F\left(L\right)$ (see appendix A for its explicit form). In the linear approximation, the normal distribution (16) implies that L is also normal, specifically

$$\begin{equation} L | L^{^{\prime}} \sim \mathcal{N}(L^{^{\prime}},\Sigma), \quad \textrm{where} \ \ \Sigma=\left[\nabla F\left(L\right)\right]^{-1}C \left[\nabla F\left(L\right)\right]^{-\mathrm{T}}.\end{equation} \tag{ 17 }$$

We deduce the probability that a LoR L is detected given that a positron is annihilated at x by marginalizing over the angles $\phi^{^{\prime}}$ and $\theta^{^{\prime}}$ of Lʹ distributed according to (14). Using that $g/R \ll 1$, we can approximate the Gaussian in $\phi^{^{\prime}}$ and $\theta^{^{\prime}}$ by $\delta(\phi^{^{\prime}}-\phi)\delta(\theta^{^{\prime}}-\theta)$ to obtain this probability as

$$\begin{equation} \frac{\sin \theta}{2\pi\left|\det\Sigma_{2}\right|^{1/2}}\exp\left(-\boldsymbol{\Delta}^{T}\Sigma_{2}^{-1}\boldsymbol{\Delta}/2\right), \end{equation} \tag{ 18 }$$

where Σ₂ is the lower-right $2\times2$ block of Σ, given explicitly in (A.10), and

$$\begin{equation} \boldsymbol{\Delta}=\left(\begin{array}{c} \Delta_{\varphi}\\ \Delta_{\theta} \end{array}\right)=\left(\begin{array}{c} a_{\varphi}-x_{\varphi}\\ a_{\theta}-x_{\theta} \end{array}\right). \end{equation} \tag{ 19 }$$

Finally, the intensity of the Poisson process in the space of LoRs and time is obtained as

$$\begin{equation} \lambda\left(\left.L\right|\boldsymbol{x},\rho\right)=\frac{\rho}{2\pi\left|\det\Sigma_{2}\right|^{1/2}}\exp\left(-\boldsymbol{\Delta}^{T}\Sigma_{2}^{-1}\boldsymbol{\Delta}/2\right) \, \chi(L \in \mathcal{D}(\boldsymbol{x}_k)) , \end{equation} \tag{ 20 }$$

with the factor $\sin \theta$ absorbed in the volume element $\mathrm{d} \mu$ in (3). The characteristic function $\chi(\cdot)$ in (20) ensures that λ is non-zero only if L belongs to the set $\mathcal{D}(\boldsymbol{x}_k)$ of LoRs through x _k that intersect the detector.

A particularly simple expression for λ is obtained under the assumption $a_{\varphi},a_{\theta}\ll R$ valid when the particle is located well away from the detector, as is the case for experiments carried out in small vessels. Equation (20) is then well approximated by setting $a_\varphi = a_\theta = 0$ to obtain

$$\begin{equation} \lambda\left(\left.L\right|\boldsymbol{x},\rho\right)=\frac{\rho}{\pi g^{2}\sin\theta}\exp\left(-\left[\frac{\Delta_{\varphi}^{2}}{g^2}+\frac{\Delta_{\theta}^{2}}{g^{2}\sin^2 \theta}\right]\right) \, \chi(L \in \mathcal{D}(\boldsymbol{x}_k)). \end{equation} \tag{ 21 }$$

Statistically, a proportion of the photons emitted from annihilation events are scattered, resulting in outlier LoR detections (see figure 1 for a visualization of such detection). We account for these outliers by adding to the Poisson process describing the generation of LoRs by K particles a component independent of the particle positions and activities. This captures the fact that scattering leads to the identification of spurious LoRs whose distribution is not directly related to the particle positions. For simplicity, we assume this distribution to be uniform over the entire space of LoRs. We therefore re-write the probability density (7) for the detected LoRs as

$$\begin{equation} P\left(\mathcal{L}|\boldsymbol{X},\varrho,\Delta t\right)=\frac{(\Delta t)^N \mathrm{e}^{-\Lambda(\boldsymbol{X},\varrho,\Delta t)}}{N!} \prod_{n=1}^N \left( \frac{\rho_0}{\sigma_0^2} + \sum_{k=1}^K \lambda\left(L_{n} |\boldsymbol{x}_k,\rho_k \right) \right), \end{equation} \tag{ 22 }$$

where

$$\begin{equation} \Lambda\left(\boldsymbol{X},\varrho,\Delta t\right)=\Delta t\left(\rho_{0}+\sum_{k=1}^{K}\int_{\mathcal{D}}\lambda\left(\left.L\right|\boldsymbol{x}_{k},\rho_{k}\right)\,\mathrm{d} \mu\right), \end{equation} \tag{ 23 }$$

in which $\mathcal{D}$ is the set of all the possible LoRs for a specific detector configuration, ρ₀ is the rate of scattered events, which is treated as independent of the activities ρ_k for simplicity, and

$$\begin{equation} \sigma_{0}^{2}=\int_{\mathcal{D}}\mathrm{d}\mu \end{equation} \tag{ 24 }$$

is the outlier variance. Note that the activities ρ_k now stand for the rate of generation of 'true' LoRs, as opposed to outlier LoR resulting from scattering events. The total activity from all particles is $\sum_{k = 0}^{K}\rho_{k}$.

The evaluation of the integrals (23) and (24) defining Λ and σ₀ requires to determine the sets $\mathcal{D}$ and $\mathcal{D}(\boldsymbol{x_k})$ of, respectively, all LoRs intersecting the detectors twice and LoRs through x _k intersecting the detectors twice. We detail in appendices B and C the parameterization of these sets and the computation of the integrals. The outlier variance is found as

$$\begin{equation} \sigma_{0}^{2}=\frac{\pi HR}{8\beta}\left(2\beta+1-\sqrt{4\beta^{2}+1}\right), \end{equation} \tag{ 25 }$$

where $\beta = R/H$. This is proportional to the detecting surface area $A_{d} = 2\pi RH$. In the limit $\beta\rightarrow0$ corresponding to a cylindrical detector that resembles a long tube, all LoRs are detected and $\sigma_{0}^2\rightarrow A_{d}/8$. In the limit $\beta\rightarrow\infty$, corresponding to a detecting surface that resembles a ring in which no true LoRs are detected, $\sigma_{0}^2\rightarrow0$ as expected. The rate Λ is approximated as

$$\begin{equation} \Lambda\left(\boldsymbol{X},\varrho,\Delta t\right)=\Delta t\left(\rho_{0}+\sum_{k=1}^{K}\rho_{k}\mathcal{G}\left(\boldsymbol{x}_{k}\right)\right), \end{equation} \tag{ 26 }$$

assuming that the distance between the particles and the detectors is much larger than g. Since $g \ll R$, this assumption holds for realistic PEPT setups. Here, $\mathcal{G}\left(\boldsymbol{x}\right)$ is a solid angle subtended by the cylindrical detector from x _k divided by 4π; it can be written as

$$\begin{equation} \mathcal{G}\left(\boldsymbol{x}\right)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\varphi\int_{\theta_{\text{min}}}^{\pi/2}\sin\theta\,\mathrm{d}\theta=\frac{1}{2\pi}\int_{0}^{2\pi}\cos\theta_{\text{min}}\left(\boldsymbol{x},\varphi\right)\,\mathrm{d}\varphi, \end{equation} \tag{ 27 }$$

where $\theta_\mathrm{min}$, given explicitly in (C.2), is the minimum angle required for the LoR to intersect the cylindrical detector twice.

According to (26), $\rho_k \mathcal{G}(\boldsymbol{x}_k)$ can be interpreted as an 'effective' activity for particle k that accounts for failure of detection of LoRs caused by detector geometry through the factor $0\leqslant \mathcal{G}(\boldsymbol{x}_k) \leqslant 1$. Figure 3 shows $\mathcal{G}$ as a function of the cylindrical coordinates (r, z) of the particle ($\mathcal{G}$ is independent of the azimuthal angle by symmetry) for three aspect ratios $\beta = R/H$. As β is increased $\mathcal{G}$ trivially decreases throughout the domain; the maximum value is always obtained at the origin $r = z = 0$. We derive asymptotic approximations for $\mathcal{G}$ at small and large β in appendix C.

Zoom In Zoom Out Reset image size

We carry out our analysis on synthetic data produced by GATE [22], a unique software for positron emission simulations. This has the benefit of providing the exact particle position as ground truth, thus enabling the exact measurement of the inference error. We consider a cylindrical detector surface with R = 200 and $H = 230~\mathrm{mm}$ ($\beta = R/H = 0.87$), composed of identical $4 \times 4$ mm detector elements. The rotating, positron-emitting particle is suspended in water, and its 'corrected' activity, that is accounting for the effect of detector deadtime which is not represented explicitly in the present work, is $\rho = 5\cdot10^4$ emissions/s. LoRs are recorded over a period of one second, providing ample data for the inference. We carry out simulations in four different conditions with varying circle radii and rotation frequencies: r = 50 mm at f = 1 and 2 Hz, and r = 100 mm at f = 0.5 and 1 Hz. This yields two distinct particle velocities $2\pi f r\approx0.3$ and 0.6 m s⁻¹, each repeated twice at varying acceleration $4\pi^2 f^2 r$.

We employ a standard Metropolis–Hastings Monte Carlo Markov Chain (MCMC) algorithm [23] to sample the posterior distribution for the particle position x and activities ρ₀ and ρ₁. The likelihood and prior are those given in (22), with λ in (20), and (10). For the current setup with $4 \times 4$ mm² detectors, the detecting error variance in (20) gives g = 2.43 mm.

The MCMC algorithm produces samples from the posterior distribution as a sequence of (statistically dependent) sets of inferred parameters. This sequence is constructed iteratively, with a new set of parameters obtained by proposing a random change to a single randomly chosen parameter, then accepting or rejecting the change depending on a likelihood ratio [21]. We use Gaussians as proposal density functions for the difference between each old and new parameter. We increase/decrease the variance of these Gaussians by a factor of two if 40 consecutive steps were accepted/rejected, which resulted in an averaged acceptance rate of 0.26. The initial values of the Gaussian proposal density are set as follows: the position variance is initially set to g²; the activity standard deviation is initially taken as 5% of the particle activity. We use 10⁵ steps to infer the position and activity at each run. We start each MCMC sampling with $\rho_{0} = \rho_{1} = \rho/2$ (with ρ the prescribed particle activity); to speed up the computation, we start from the exact particle position thus not requiring a burn-in period for the MCMC.

For each pair of values of the radius r and frequency f of the particle trajectory, we infer the position and activity of the particle at ten equally spaced times t_i using the LoRs collected in the time interval $t_{i}-\Delta t/2\lt t_{i}\lt t_{i}+\Delta t/2$. The statistics of each set of 10 MCMC samplings makes it possible to characterize the variability of the inference process that arises from the randomness of the positron emission. Our main focus is on the dependence of the accuracy of the inference on the choice of the observation interval $\Delta t$. We therefore run the entire MCMC computation for a range of values of $\Delta t$ and compare the results.

We first show the results of a single MCMC sampling, carried out for a particle rotating with frequency $f = 1~\mathrm{Hz}$ at a radius $r = 50~\mathrm{mm}$, resulting in the tangential velocity $v = 2\pi f r\approx0.3~\mathrm{m~s}^{-1}$. The time interval is chosen as $\Delta t = 10~\mathrm{ms}$ which yields approximately 250 detected LoRs.

In figure 4 we show two-dimensional histograms of the posterior pdf for the particle position in the $(x,y,0)$ horizontal (a) and $(x,-3.56,z)$ vertical (b) planes. (the particle true position at the midst of the time interval is $\boldsymbol{x} = \left(49.87,-3.56,0\right)$) The results illustrate the localization of the posterior defined by (22), with nearly the entire pdf supported in a $2\times 2 ~ \mathrm{mm}^2$ area in both the horizontal and vertical. To quantify the uncertainty in the inference of the particle position more systematically, we estimate a covariance matrix for the particle position from the MCMC samples. We use its three eigenvalues, $\mathscr{L}_i$ say, to compute the metric

$$\begin{equation} s=\mathcal{R}_{3}\left(\prod_{i=1}^{3}\sqrt{\mathscr{L}_{i}}\right)^{1/3}, \end{equation} \tag{ 28 }$$

where $\mathcal{R}_{3} = 2.79$, such that s is the radius of a sphere that marks the 95% confidence level for the Gaussian distribution with the estimated covariance. ($\mathcal{R}_{3}$ is the three-dimensional equivalent of the more famous $\mathcal{R}_{1} = 1.96$ that gives the 95% confidence interval for a one-dimensional Gaussian distribution.) The pdf in figure 4 has s = 0.79 mm, five times smaller than the detector side of 4 mm.

Zoom In Zoom Out Reset image size

Next, we run multiple MCMCs with increasing $\Delta t$ in order to assess the trade-off between (i) the reduction in uncertainty associated with the larger number of LoRs, and (ii) the increase in uncertainty caused by the particle displacement during $\Delta t$. Indeed, in the absence of motion, we expect s to decrease as $N^{-1/2} \approx (\rho\Delta t)^{-1/2}$, but the particle displacement introduces an error in the model underpinning (22) that increases with $\Delta t$. One may wonder why the activity ρ is not increased to alleviate the uncertainty in the moving particle position. Increasing ρ provides more data over the same $\Delta t$ and hence monotonically decreases the uncertainty. In fact, ρ may be sufficiently increased such that the particles are quasi-static over $\Delta t$, in which case $s\propto \rho^{-1/2}$. However, in PEPT, as opposed to PET, activity cannot be increased by simply introducing more radioactive material with the same activity density. Here each particle is radiolabeled separately, and the requirement that particles remain small so that they are easily suspended in fluid makes it challenging to increase each particle activity. We therefore analyze particles with fixed activity, rotating with the two tangential velocities 0.3 and 0.6 m s⁻¹ obtained by changing either the frequency or radius of the circular trajectory. In this way, we can examine the effect of varying the acceleration while maintaining a constant velocity. The results are shown in figures 5(a) for v = 0.3 m s⁻¹ and 5(b) for v = 0.6 m s⁻¹. The blue curves, with the left vertical axis, show the metric s, averaged over the ten positions of the particle, as a function of $\Delta t$, with the error bars indicating the corresponding standard deviation. The red curves, with the right vertical axis, show the absolute error $\left|\boldsymbol{x}-\boldsymbol{x}_*\right|$ between the average inferred position x and the particle actual position $\boldsymbol{x}_*$. The solid and dashed curves are calculations for particles rotating with r = 50 and 100, respectively, corresponding to a twofold decrease in particle acceleration. The scales in both the left and right vertical axes in figure 5 are equal for ease of comparison.

Zoom In Zoom Out Reset image size

All the uncertainty curves (blue, left vertical axis) exhibit minima, reflecting the trade-off mentioned above. As could be expected, s is smaller for v = 0.3 m s⁻¹ than for v = 0.6 m s⁻¹ as a result of smaller particle displacements with the same rate of LoR generation, and for the same velocity, s is smaller for the smaller acceleration. The curves are asymmetric about the minima, with a sharp decrease, and a much milder increase (and even some further decrease at large $\Delta t$ for v = 0.6 m s⁻¹ and r = 100 mm (solid curve in figure 5(b))). The error (red, right vertical axis), in contrast, increases unambiguously for large $\Delta t$, reflecting the expected model's inability to accurately infer the instantaneous position when the particle moves significantly over $\Delta t$. To provide reference to our results, we also analyzed the data with the Birmingham algorithm [1], the most commonly used method for PEPT inversion [14], and the expectation-maximization algorithm [9]. The error $\left|\boldsymbol{x}-\boldsymbol{x}_*\right|$ achieved by both methods (not shown in the figures) is roughly 1 mm—typically larger than our error at small $\Delta t$ but smaller at large $\Delta t$. Other inversion algorithms may be tested and compared against our results; however, our primary aim is to quantify the uncertainty in PEPT rather than developing an alternative algorithm. Accordingly, no attempts were made to optimize our method by, for instance, testing various priors.

It might seem paradoxical that the uncertainty does not increase substantially with increasing $\Delta t$ while the error does. This is the result of our weakly constrained model of scattering: LoRs that are inconsistent with emission by a fixed particle are attributed to the scattered component of the likelihood, whether they arise because of scattering or because of particle motion. This explanation is confirmed in figure 6 which shows the inferred relative activity $\psi = \rho_1/\left(\rho_0+\rho_1\right)$ of unscattered LoRs (blue curves, left axis). Although, the correct, physical activity should be independent of $\Delta t$, the inferred value decreases, as a result of the increasing misattribution of LoRs to the scattered component. One way of remedying this artifact is to constrain the ratio $\rho_1/\rho_0$ or, equivalently, ψ to a physically sensible range using the prior. Another way is to account for the particle motion in the likelihood. We adopt the latter approach next by implementing the Taylor-expansion-based first-order inference described in section 2.2.

Zoom In Zoom Out Reset image size

The first-order inference relies on the form (11) of Bayes' formula, with the time-dependent likelihood (4) and a linear approximation to the particle trajectory, so both the particle position and velocity in the middle of the observation interval need to be inferred. For simplicity, we approximate $\boldsymbol{x}(t) \approx \boldsymbol{x}(t_0) = \mathrm{const.}$ in the factor Λ in (5) and use expression (8) that depends on positions only. The dependence on velocities is then through the rates $\lambda(L|\boldsymbol{x}_k(t),\rho_k)$.

We carry out MCMC sampling for the first-order (M = 1) inference in the case r = 100 mm and v = 0.6 m s⁻¹ and compare it with the results of the position-only inference (M = 0) of the previous section. For this comparison we use the same number of MCMC steps for M = 1 as for M = 0, even though the number of parameters to estimate increases from 5 to 8. The results are shown in figure 7(a) which shows the uncertainty metric s (blue, left axis) and the position error (red, right axis), with the same conventions as in figure 5.

Zoom In Zoom Out Reset image size

Both the M = 1 uncertainty and error curves lay below their counterparts for M = 0. This demonstrates that inferring the instantaneous velocity increases the measurement accuracy. Furthermore, the variance of the set of ten experiments, indicated by the error bars, is also smaller. The difference between the errors for M = 0 and 1 is marginal at $\Delta t \lesssim 20$ ms, where the particle moves less than $v\Delta t \approx 12$ mm over the time interval. At larger $\Delta t$, however, the M = 0 error curve rises at a much steeper rate compared with the M = 1 one. The increased particle movement is then tackled better by inferring the particle trajectory over $\Delta t$ as linear in t rather than fixed. Note that we carried out the M = 1 inference for $\Delta t$ up to 200 ms (not shown in figure 7(a)) to verify that s reaches a minimum value, found to be $\Delta t\approx150$ ms.

A benefit of M = 1 inference over M = 0, in addition to an increase accuracy, is that it estimates directly the particle velocity, often the main quantity of interest for fluid dynamical applications, which otherwise needs to be deduced from the positions, e.g. using finite differences. We show the uncertainty and error on the inferred particle velocity in figure 7(b). The blue curve (left axis) shows the mean (over ten particle positions, as before) uncertainty metric s_v defined as in (28) but using the eigenvalues of the velocity covariance matrix, with error bars indicating the standard deviation. The red curve shows the mean error $|\boldsymbol{v}-\boldsymbol{v}_*|$ between the inferred particle velocity v and the exact velocity of the circular motion $\boldsymbol{v}_*$. The qualitative behavior of the uncertainty and error on the inferred velocity is roughly similar to that corresponding to the position, though the minimum error appears for a larger $\Delta t$. The decrease in inferred uncertainty as $\Delta t$ increases is related to the misattribution of LoRs to the scattered components, as discussed in section 4.3.