An analytic, moment-based method to estimate orthopositronium lifetimes in positron annihilation lifetime spectroscopy measurements

Abstract

The presence of tumor hypoxia is known to correlate with poor patient prognosis. Measurement of tissue oxygen concentration can be challenging, but recent advancements using positron annihilation lifetime spectroscopy (PALS) in three-dimensional positron emission tomography (PET) scans have shown promise for hypoxia detection. In this work, a novel method for estimating the orthopositronium lifetime in PALS is presented. This method is analytical and uses moments of the time-difference histogram from photon arrival times. For sufficient statistical power, the method produces monotonic, stable estimates. For cases with a lower number of photon counts, the method was characterized and solutions are presented to correct for bias and estimation variability.

I Introduction

Positronium lifetime imaging (PLI) is a novel augmentation for positron emission tomography (PET) that may allow PET scans to extract tissue oxygenation information, including hypoxic locations, in addition to the specific biochemical properties the employed PET tracer is responsive to. Recently, Moskal et al. has identified the positronium lifetime as a possible metric for hypoxia [1]. The same group has also designed and built a novel PET scanner for generating positronium lifetime images [2]. Other recent work includes that of Shibuya et al. in which the correlation of oxygen concentration with positronium lifetimes was rigorously established [3].

PLI makes use of measurements of the lifetime of positronium (Ps), which are short-lived bound states between a positron emitted from nuclear decay and an electron from the environment. In clinical PET circumstances, approximately 40% of the released positrons lead to this electron-positron bound state [3]. The other 60% annihilate with electrons without forming positronium. These produce direct annihilations (DA), and predominantly produce two coincident $511\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ photons with a lifetime of about $388\text{\,}\mathrm{p}\mathrm{s}$ [2].

There are two types of positronium, orthopositronium (o-Ps) where the electron and positron spins are parallel, and parapositronium (p-Ps) where the spins are anti-parallel. Both are formed within $5\text{\,}\mathrm{p}\mathrm{s}$ of the release of the positron after it has thermalized with the environment. In tissue, approximately 75% of positronium is o-Ps and 25% is p-Ps. Both o-Ps and p-Ps are inherently unstable and eventually their constituent electron and positron annihilate with each other without interacting with their environment. In such self annihilations, p-Ps produces two $511\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ photons and has a short lifetime of approximately $125\text{\,}\mathrm{p}\mathrm{s}$. On the other hand, o-Ps produces three coincidence photons and has a long lifetime of approximately $142\text{\,}\mathrm{n}\mathrm{s}$ [4].

As a result, o-Ps has the opportunity to interact with electrons that are present in its surroundings and to decay before self-annihilation, thereby shortening the lifetime to 1-$15\text{\,}\mathrm{n}\mathrm{s}$ [4]. These interactions include a spin-exchange with an unpaired electron that converts o-Ps to p-Ps where it quickly decays to produce two $511\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ photons, and a pick-off event where the positron annihilates with a local electron (also most likely to produce two $511\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ photons). The most dominant of these interactions is the spin-exchange. In tissue, the o-Ps lifetime is between 1-$3\text{\,}\mathrm{n}\mathrm{s}$ and importantly its value depends on local properties that affect the strength of the spin-exchange interaction [3]. For example, molecular oxygen is paramagnetic and can readily cause o-Ps decay, since o-Ps readily decays by interacting with unpaired electrons.

If we therefore choose an isotope that possesses a nontrivial decay mode in which a prompt gamma and a positron are released at essentially the same time, we may measure the arrival time difference between the prompt gamma and coincident $511\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ photons. A metric for local tissue oxygenation may then be derived. Specifically, the histogram $S(\Delta t)$ of these differences in arrival times $\Delta t$ can be modeled as the sum of multiple convolutions of exponential distributions and Gaussian distributions [2]. We consider here lifetime components from o-Ps, p-Ps, and DA as they have sufficiently distinct lifetimes and intensities to be modeled separately. Each is represented by one exponential distribution that models the probability distribution for the decay, convolved with a Gaussian distribution that models the statistical uncertainty in time measurement. Therefore, the full model for $S(\Delta t)$ can be written as

$$S(\Delta t)=b+\sum^{3}_{i=1}I_{i}\ \textrm{EXP}({\Delta t};{\lambda_{i}})\ast\mathcal{N}({\Delta t-t_{0}};{\sigma_{i}}),$$ (1)

where $\ast$ denotes convolution, $I_{i}$, $\lambda_{i},$ and $\sigma_{i}$ are the intensity, decay-rate constant, and standard deviation of the time-measurement uncertainty associated with component $i$ for $i=1,2,3$, $b\geq 0$ accounts for the presence of background events, EXP is defined as

$$\textrm{EXP}({t};{\lambda})=\lambda\exp\{-\lambda t\}\cdot u(t),$$ (2)

in which $u(t)$ is the unit step function defined by $u(t)=1$ for $t\geq 0$ and $u(t)=0$ for $t<0$, and

$$\mathcal{N}({t};{\sigma})=\left(\sqrt{2\pi}\sigma\right)^{-1}\,\exp\{-t^{2}/2\sigma^{2}\}.$$ (3)

The lifetime $\tau_{i}$ of component $i$ is related to $\lambda_{i}$ by $\tau_{i}=1/\lambda_{i}$. In Equation 1, $t_{0}$ is introduced to allow for an offset in time measurement. We wish to estimate $\tau$ associated with o-Ps from a measurement of $S(\Delta t)$ for providing a metric for certain tissue property. Typically, it is reasonable to take $\sigma_{i}=\sigma,\,\forall i$ and $t_{0}=0$. Without loss of generality, we can take the components $i=1,2,3$ to be those from the p-Ps decay, DA, and o-Ps decay, respectively, and hence $\lambda_{1}>\lambda_{2}>\lambda_{3}.$

While PLI is emerging for the purpose of producing full three-dimensional images of the o-Ps lifetime, this work is concerned with only the single-dimensional case represented by Equation 1. For distinction, we therefore refer to this method as one for positron annihilation lifetime spectroscopy (PALS) measurements. In the discussion we will comment on using this method for imaging in general, and will consider it in-depth in a separate paper.

Presented here is a novel analytic method for estimating the o-Ps lifetime from $S(\Delta t)$. Current methods for this task include fitting a single exponential distribution to the time-difference histogram [3], and fitting a reduced full model given by Equation 1 by assuming known values for certain parameters to the histogram [2]. Additionally, Shibuya et al. has proposed an inverse Laplace transform method to distinguish between positronium lifetimes while merging voxels for better statistics [5]. Their method is able to discern similar lifetime values, but still employs curve-fitting. In comparison with these curve-fitting methods, the new analytic method in this work is more computationally efficient, which is an important consideration for future application of the method for PLI in which lifetimes need to be obtained for a large number of voxels in a clinical setting. In addition, under mild and realistic conditions the analytic method is not sensitive to the unknown lifetimes for DA and p-Ps, nor to the time-measurement uncertainty $\sigma$. Generally, this is the not case with curve-fitting methods.

II Materials and Methods

II.1 Derivation of the method and justifications

The $n$th moment of a function $f(t)$ is defined by

$$\mu_{n}\{f(t)\}=\int_{-\infty}^{\infty}t^{n}f(t)dt$$ (4)

if the integral exists. The proposed analytic method is based on the following observations.

1. The $n$th moment of $\textrm{EXP}({t};{\lambda})\}$ exists for all $\lambda>0$ and is given by

   $$\mu_{n}\{\textrm{EXP}({t};{\lambda})\}=\frac{n!}{\lambda^{n}}.$$ (5)

2. Let $g(t;\lambda,\sigma)=\textrm{EXP}({t};{\lambda})\ast\mathcal{N}({t};{\sigma}).$ For sufficiently large $\lambda$ (with respect to $\sigma$) and $n$, it can be shown that, for any $s>0$,

   $$\mu_{n}\{e^{-st}g(t;\lambda,\sigma)\}\approx c\cdot\mu_{n}\{\textrm{EXP}({t};{s+\lambda})\},$$ (6)

   where $c=\lambda/(s+\lambda)e^{\sigma^{2}\lambda^{2}/2}$ is independent of $n$.

3. Applying the above to Equation 1 with $\sigma_{i}=\sigma$ and $t_{0}=0$, we get

   $$\mu_{n}\{e^{-s\Delta t}\tilde{S}(\Delta t)\}\approx c\cdot n!\cdot\sum_{i=1}^{3}\frac{I_{i}}{(s+\lambda_{i})^{n}},$$ (7)

   where $\tilde{S}(\Delta t)=(S(\Delta t)-b)\cdot u(\Delta t)$.

4. Since $\lambda_{3}$ is assumed to be smaller than $\lambda_{1}$ and $\lambda_{2}$, the $i=3$ term dominates the sum in Equation 7 for large values of $n$, yielding

   $$R(n+1,s)\equiv\frac{\mu_{n+1}\{e^{-st}\tilde{S}(\Delta t)\}}{\mu_{n}\{e^{-st}\tilde{S}(\Delta t)\}}\approx\frac{n+1}{s+\lambda_{3}}.$$ (8)

   Therefore, with an adequately large $n$ we have

   $$\tau_{3}\approx\frac{R(n+1,s)}{n+1-sR(n+1,s)+\delta},$$ (9)

   where the small positive value $\delta$ has been added to the denominator to control observed estimation variability for small values of $\tau_{3}$.

The plot in Figure 1 demonstrates Equation 6 numerically for $s=1$ and several selected values for $n$ and $\lambda$. Derivations of Equation 5 and 6 are given in the Appendix.

Figure 1: The scale factor $c$ in Equation 6, from which the ratio ${\mu_{n}\{e^{-st}g(t;\lambda,\sigma)\}}/{\mu_{n}\{\textrm{EXP}({t};{s+\lambda})\}}$ has been subtracted. Their difference has a maximum of $2.01\text{\times}{10}^{-03}\text{\,}$, which decreases as $n$ increases until $n=12$ where the mean difference becomes $1.6\text{\times}{10}^{-16}\text{\,}$.

For the constant values we used in Equation 1, we reference Moskal et al. [2]. These values are $\tau_{1}=1/\lambda_{1}\approx$ $0.125\text{\,}\mathrm{n}\mathrm{s}$ and $I_{1}\approx 0.1$ for p-Ps, $\tau_{2}=1/\lambda_{2}\approx$ $0.388\text{\,}\mathrm{n}\mathrm{s}$ and $I_{2}\approx 0.6$ for DA, and $\tau_{3}=1/\lambda_{3}\approx$ 1-$10\text{\,}\mathrm{n}\mathrm{s}$ and $I_{3}\approx 0.3$ for o-Ps. The coincidence resolving time (CRT) was chosen to be $600\text{\,}\mathrm{p}\mathrm{s}$ full width at half maximum (FWHM) from current state-of-the-art time-of-flight (TOF) PET systems [6, 2], yielding $\sigma\lesssim$$220\text{\,}\mathrm{p}\mathrm{s}$. Therefore, for PALS and PLI measurements, the assumptions leading to Equation 9 can be verified.

In theory, $\tilde{S}(\Delta t)$ decreases exponentially to zero as $\Delta t$ increases, which allows $\mu_{n}\{\tilde{S}(\Delta t)\}$ to be well defined. In reality, noise in $\tilde{S}(\Delta t)$ does not necessarily decrease with $\Delta t$ and hence will contribute a substantial statistical error in $\mu_{n}\{\tilde{S}(\Delta t)\}$, especially for large $n$. This instability can be considerably reduced by instead using $\mu_{n}\{e^{-s\Delta t}\tilde{S}(\Delta t)\}$, $s>0$ as the term $e^{-s\Delta t}$ attenuates the contribution from data at large $\Delta t$. Although a large $s$ is favored for alleviating the effects of noise, it also diminishes the differences among $s+\lambda_{i}$ and thereby requires a large $n$ for Equation 8 to hold true. However, calculation of higher-order moments is more susceptible to noise so the value of $s$ needs to be chosen with care. In this work, in consideration of the numerical values of $\lambda_{i}$ and $I_{i}$, we choose $s=1$, which is empirically justified based on the results to be reported below. Future work will consider $s$ more extensively.

From Equation 4, $\mu_{n}\{f(t)\}$ can be regarded as a filter that removes an increasingly wider range of small $t$ data as $n$ increases. By inspection of Figure 2 it can be stated that by choosing an $n$ which makes Equation 8 valid, a “soft” low cutoff has been introduced. This avoids using data where DA and p-Ps are dominant.

II.2 Simulated data

The proposed method is evaluated by using computer generated simulation data. All computation codes were implemented in Python version 3.11, using specified values for $b$, $I_{i}$, $\lambda_{i}$, and $\sigma$, and the desired total number of events for the histogram $N$. The simulation program first computed $S(\Delta t_{i})$ according to Equation 1 at discrete times in [$-5\text{\,}\mathrm{n}\mathrm{s}$, $25\text{\,}\mathrm{n}\mathrm{s}$] and at a regular spacing of $60\text{\,}\mathrm{p}\mathrm{s}$. Next, it scaled $S(\Delta t_{i})$ to yield $\sum_{i}S(\Delta t_{i})=N$. Then, the scaled $S(\Delta t_{i})$ was replaced by an integer drawn by a Poisson random number generator (from the numpy.random module) whose mean equals the scaled $S(\Delta t_{i})$. Figure 2 shows an example of the generated histogram where the parameters were so chosen that it was similar to the measured histogram reported by Moskal et al. [2]. Each simulated histogram, except where noted, contained $N=$2\text{\times}{10}^{5}\text{\,}$$ total events.

Figure 2: A simulated histogram of $\Delta t$ having $2\text{\times}{10}^{5}\text{\,}$ events and using a $60\text{\,}\mathrm{p}\mathrm{s}$ bin size. It was generated by using the following parameters: $\lambda_{1}I_{1}$, $\lambda_{2}I_{2}$, $\lambda_{3}I_{3}$ = $0.078$, $0.388$, $0.165$, $\tau_{1}$, $\tau_{2}$, $\tau_{3}=$0.125\text{\,}\mathrm{n}\mathrm{s}$$, $0.388\text{\,}\mathrm{n}\mathrm{s}$, $2.5\text{\,}\mathrm{n}\mathrm{s}$, $\sigma=$0.169\text{\,}\mathrm{n}\mathrm{s}$$, $b=10$, and $t_{0}=$9\text{\,}\mathrm{n}\mathrm{s}$$. $\lambda_{1},\lambda_{2},\lambda_{3}=$ $8\text{\,}\mathrm{n}\mathrm{s}^{-1}$, $2.58\text{\,}\mathrm{n}\mathrm{s}^{-1}$, $0.4\text{\,}\mathrm{n}\mathrm{s}^{-1}$. The individual contributions from the p-Ps, DA, and o-Ps components are shown in yellow, green, and magenta, respectively.

II.3 Implementation and numerical studies

The background term $b$ was estimated using the average of the histogram in the $\Delta t<$0\text{\,}\mathrm{n}\mathrm{s}$$ region. The estimate was then subtracted from each histogram data points and $\lambda_{3}$ was calculated according to Equation 8 and Equation 9. For computing moments, Equation 4 was approximated by naive discrete summation: $\mu_{n}\{f(t)\}\approx\delta t\sum_{i}t_{i}^{n}\tilde{S}(t_{i})$ where $t_{i}$ are the time points where $f(t)$ is available and $\delta t$ is the spacing of $t_{i}$.

The proposed method was evaluated for accuracy and precision against a number of parameters, including 1) the order of moment $n$ used, 2) the upper truncation of the histogram data, 3) the number of counts $N$, and 4) the background level $b$. At present, the range of in vivo o-Ps lifetime values has not been precisely established. However, in a recent paper Moskal et al. observed that cardiac myxoma and adipose tissue had mean lifetimes of $1.912\text{\,}\mathrm{n}\mathrm{s}$ and $2.613\text{\,}\mathrm{n}\mathrm{s}$, respectively [2]. Therefore, we performed evaluations for o-Ps lifetimes of 1.0, 1.5, 2.0, 2.5, and 3.0 ns to cover the likely in vivo lifetime range. On the other hand, since they are insensitive to the local environment [7], the reported mean values of $388\text{\,}\mathrm{p}\mathrm{s}$ and $142\text{\,}\mathrm{p}\mathrm{s}$ are used for DA and p-Ps lifetimes, respectively. For $I_{1}$, $I_{2}$ and $I_{3}$, the values used were based on quantitatively matching the simulated and measured data.

III Results

The results are presented by Figure 3 through Figure 10. Each figure contains individual points which are estimates of $\tau_{3}$, denoted by $\hat{\tau}_{3}$ below, derived from our method. Each data point is the mean of the results obtained from $1\text{\times}{10}^{4}\text{\,}$ histograms simulated by using the same parameters ($1\text{\times}{10}^{5}\text{\,}$ in the case of Figure 3), and the shaded regions in the plots give the $\pm 1$ standard deviations (std) about the means. The horizontal lines, when present, indicate the true o-Ps lifetimes that are used to produce simulation data.

III.1 Lifetime estimate versus order of moment $n$

Figure 3 shows the estimated o-Ps lifetime when the order of moment $n$ employed by Equation 9 is varied. Four general trends can be observed. First, all curves show a plateau where the estimated lifetime has essentially zero bias. This plateau occurs between $n\approx 5$ and $n\approx 16$ depending on the o-Ps lifetime. The standard deviations of the estimates are also sufficiently small to allow for discrimination of all the lifetimes examined. Second, the standard deviation increases with $n$, which is consistent with the observation made above that higher-order moments are more sensitive to data noise. Third, as $\tau_{3}$ increases the plateau occurs at a lower $n$. This is because the differences between $s+\lambda_{3}$ and $s+\lambda_{i}$, $i=1,2$ increases, allowing the $i=3$ term to dominate the sum in Equation 7 at a smaller $n$. Fourth, all curves decrease toward zero as the order increases. This is because higher-order moments are increasingly more contributed by data at larger $\Delta t$ while data is simulated only for $-$5\text{\,}\mathrm{n}\mathrm{s}$\leq\Delta t\leq$25\text{\,}\mathrm{n}\mathrm{s}$$. In practice, the measured histogram is necessarily truncated.

Figure 3: Estimated o-Ps lifetime as a function of the order of moment $n$. The shaded areas for each curve represent the $\pm 1$ standard deviation on each data point. The true $\tau_{3}$ values shown in the legend are those used for producing simulation data.

It is also noted that all curves converge at small $n$. This reflects the situation that with a small $n$ the $i=1,2$ terms in fact dominate the $i=3$ term in Equation 7. As a result, Equation 9 yields some fixed value since $I_{1,2}$, and $\tau_{1,2}$ are constants. Based on the plot, we propose to that in general $n\lesssim 5$ should not be used for PALS o-Ps estimations. Table III.1 summarizes the bias-minimizing order of moment in the plateau region and the statistics of the estimates obtained when using this order. As shown, the largest std/mean ratio is only 0.654%, obtained for $\tau_{3}=$1.0\text{\,}\mathrm{n}\mathrm{s}$$.

Table 1: Bias-minimizing estimates of Figure 3, shown with the moment $n$ which minimized the bias. We again emphasize that this moment-based method requires $n$ to be sufficiently large.

$\tau_{3}$ (ns)	order $n$	$\hat{\tau}_{3}$
		mean (ns)	std (ns)	std/mean	% error
1.0	11	0.99	0.094	0.095	0.654
1.5	11	1.50	0.056	0.037	0.096
2.0	12	2.00	0.051	0.025	0.008
2.5	10	2.50	0.041	0.016	0.001
3.0	10	3.00	0.037	0.012	0.009

Figure 4 shows the result obtained when the $\tau_{3}$ values measured for cardiac myxoma and adipose tissue by Moskal et al. [2] were used to produce simulation data. In this case, the bias-minimizing orders of the moments were found to be $n=11$ for $\tau_{3}=$1.912\text{\,}\mathrm{n}\mathrm{s}$$ (myxoma tissue), which yielded $\hat{\tau}_{3}=1.91\pm 0.05$ ns, and $n=11$ for $\tau_{3}=$2.613\text{\,}\mathrm{n}\mathrm{s}$$ (adipose tissue), which yielded $\hat{\tau}_{3}=2.61\pm 0.04$ ns. The means of these estimates are within $0.02\%$ and $0.005\%$ of their respective true values.

Figure 4: Similar to Figure 3, this shows two practical cases for o-Ps lifetimes which have been measured recently in vivo [2].

III.2 Lifetime estimate versus number of events $N$

The results discussed were obtained by using simulated data with $N=$2\text{\times}{10}^{5}\text{\,}$$. Generally, a histogram derived from a larger number of events has better statistics and can lead to better estimates for $\lambda_{3}$. Figure 5 shows the results for $N$ ranging from $1\text{\times}{10}^{3}\text{\,}$ and $1\text{\times}{10}^{8}\text{\,}$. The order of moment was fixed to $n=11$. A tabulated summary of these results are shown in Table III.2.

Figure 5: Estimated o-Ps lifetimes as a function of the number of events $N$ in the histogram. The shaded areas again indicate $\pm 1$ standard deviation about the mean.

Table 2: Selected estimates from Figure 5.

$\tau_{3}$ (ns)	$\hat{\tau}_{3}$ (ns)
	$N=${10}^{3}\text{\,}$$	$N=${10}^{4}\text{\,}$$	$N=${10}^{5}\text{\,}$$	$N=${10}^{6}\text{\,}$$
1.0	$0.6\pm 0.6$	$0.9\pm 0.4$	$1.0\pm 0.1$	$0.997\pm 0.004$
1.5	$2\pm 16$	$1.4\pm 0.3$	$1.49\pm 0.08$	$1.498\pm 0.002$
2.0	$2.0\pm 0.8$	$2.0\pm 0.2$	$2.00\pm 0.06$	$2.000\pm 0.002$
2.5	$2.5\pm 0.6$	$2.5\pm 0.2$	$2.50\pm 0.06$	$2.500\pm 0.001$
3.0	$3.0\pm 0.6$	$3.0\pm 0.1$	$3.00\pm 0.06$	$3.000\pm 0.001$

It can be seen that, as expected, for all the $\tau_{3}$ examined the standard deviations of the estimates decrease continually as $N$ increases. The estimates asymptotically approach the true values as $N$ increases and the statistics are seen to improve. The occurrences of some exceptionally large standard deviations reflect the instability of the ratio estimate given by Equation 9 when the denominator is erroneously small with respect to the numerator due to data noise. Figure 6 shows a simulated histogram for $N=$1.5\text{\times}{10}^{3}\text{\,}$$ and $\tau_{3}=$1.5\text{\,}\mathrm{n}\mathrm{s}$$ when one such case occurs. It can be seen that the $\Delta t>4$ region of the histogram where o-Ps is assumed to dominate is sparsely populated, leading to large and inconsistent errors in the numerator and denominator of Equation 9. Controlling this variability will be commented on in Section IV.3.

Figure 6: A simulated histogram for $\tau_{3}=$1.5\text{\,}\mathrm{n}\mathrm{s}$$ and $N=1,500$, including the noise-free contributions of DA, p-Ps, and o-Ps.

III.3 Lifetime estimate versus cutoff and background

In an attempt to further alleviate the deleterious effects of noise, we also examined removing data above certain $\Delta t$. Figures 7 shows that, for the case of $n=11$, as the upper truncation threshold $\Delta t_{up}$ is decreased, the standard deviation decreases, which is particularly evident for small $\tau_{3}$, whereas the bias increases in the form of underestimation. However, as shown in Figure 8, the standard deviations relative to the distances between means increases as $\Delta t_{up}$ is lowered, at least for $\Delta t_{up}$ above $\sim$10\text{\,}\mathrm{n}\mathrm{s}$$. Therefore, at least for discrimination tasks, it is beneficial to apply truncation even though it introduces bias.

Figure 7: Estimated o-Ps lifetime as a function of the upper truncation threshold $\Delta t_{up}$. Note that here $\hat{\tau}_{3}/\tau_{3}$ is plotted. The horizontal line indicates the perfect estimate given by $\hat{\tau}_{3}/\tau_{3}=1$.

Figure 8: The ratio of the standard deviation to the mean of the curves from Figure 7. The ratios are seen to decrease as the upper truncation decreases. The standard deviations about the mean values are too small to be visible.

The presence of a nonzero background $b$ increases data noise in two manners. First, for a given number of total events the number of events contributing to the signal is lowered. Second, although $b$ may be estimated and subtracted from $S(\Delta t)$, the noise associated with it remains in the data. Due to its Poisson nature, the variance of the associated noise is proportional to $b$. Figure 9 shows the relationship between the estimate $\hat{\tau}_{3}$ and the true $\tau_{3}$ for three background levels with $n=11$, $N=$2\text{\times}{10}^{5}\text{\,}$$, and no upper truncation. The background is modeled as Poisson noise with means of $b=0,10,$ and $20$ counts.

Figure 9: Estimated o-Ps lifetime versus the true value for three background levels $b$ with $n=11$, $N=$2\text{\times}{10}^{5}\text{\,}$$ and no upper truncation. The region close to $\tau_{3}=$0.5\text{\,}\mathrm{n}\mathrm{s}$$ may not be monotonic for some cases due to the curvature.

III.4 Monotonicity

We have observed that for discrimination tasks it is beneficial to allow for for some bias if the std/mean ratio of the estimate decreases. For quantitative tasks, as long as $\hat{\tau}_{3}$ is related to $\tau_{3}$ monotonically in the mean, this observation remains true because we can correct for the bias in $\hat{\tau}_{3}$ by using, for example, a predetermined calibration curve. Such monotonicity is observed in both Figure 7 and Figure 9. Since curves in Figure 3 do not cross one another this monotonicity is also true when varying $n$. For $b=0,10,$ and $20$ mean counts, $\hat{\tau}_{3}$ remains monotonic with $\tau_{3}$ in the mean and is approximately linear for $\tau_{3}\gtrsim$1.5\text{\,}\mathrm{n}\mathrm{s}$$. Below $\sim$0.5\text{\,}\mathrm{n}\mathrm{s}$$, monotonicity may not hold for $n=11$, however this is well below the range for currently-known biological values of $\tau_{3}$.

IV Discussion

IV.1 Effect of the number of events

It was observed that our moment-based method is significantly affected by the presence of Poisson noise, which comes with reduced event numbers. However, in the mean this method is still monotonic with the true value. A clinical PET system would produce a single histogram for each voxel (even though this is not necessary for estimating $\tau_{3}$ using our method). A typical scan may consist of approximately $185\text{\,}\mathrm{M}\mathrm{B}\mathrm{q}$ $($5\text{\,}\mathrm{m}\mathrm{C}\mathrm{i}$)$ of injected activity and a $64\text{\,}\mathrm{m}\mathrm{m}^{3}$ voxel size ($4\text{\,}\mathrm{m}\mathrm{m}$ side length). Assuming the human body to be mainly composed of water with a total mass of $80\text{\,}\mathrm{k}\mathrm{g}$, it would have approximately $8\text{\times}{10}^{4}\text{\,}\mathrm{c}\mathrm{m}^{3}$ of internal volume. The activity in this volume would then be $2.3\text{\times}{10}^{3}\text{\,}\mathrm{B}\mathrm{q}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$ and we might expect $148\text{\,}\mathrm{B}\mathrm{q}\mathrm{/}\mathrm{v}\mathrm{o}\mathrm{x}\mathrm{e}\mathrm{l}$ if there is a uniform activity distribution. Assuming a per-voxel sensitivity of $1\%$, 1.48 counts would be collected each second. To obtain ${10}^{3}\text{\,}$ counts for each voxel, we would need to collect data for approximately 11.3 minutes. For reference, Moskal et al. reported human brain PLI data, where the scan time was 10 minutes (after a standard radiotracer distribution waiting period), and the number of counts collected were 342 for the healthy brain tissue, 547 for the tumor, and 1119 for the salivary glands [8].

Additionally, above we have assumed a uniform radiotracer distribution. When using tumor-specific radiotracers, the radioactivity can be concentrated to tumor regions which will improve the statistics.

IV.2 Short o-Ps lifetimes in vivo

Thus far, the shortest published o-Ps lifetime value using a phantom was $1.8239\text{\,}\mathrm{n}\mathrm{s}$ [3]. This was measured in fully $\mathrm{O}_{2}$ saturated water. In a more artificial environment, Stepanov et al. [9] bubbled oxygen gas through water and measured an o-Ps lifetime of $1.746\text{\,}\mathrm{n}\mathrm{s}$. With argon gas bubbled through water, however, a lifetime of $1.833\text{\,}\mathrm{n}\mathrm{s}$ was measured. This is evidence that even extremely well-oxygenated in vivo environments will not have o-Ps lifetimes of less than $1.746\text{\,}\mathrm{n}\mathrm{s}$. In a recent review, Moskal and Stȩpień concluded that the mean biological o-Ps lifetime would be approximately $2\text{\,}\mathrm{n}\mathrm{s}$ [4].

This moment-based method has been shown to produce viable results for lifetimes above approximately $1\text{\,}\mathrm{n}\mathrm{s}$, depending on noise, moment, and truncation. For lifetimes on that are on the order of $1\text{\,}\mathrm{n}\mathrm{s}$, standard deviations for $n>15$ (seen in Figure 3) grow rapidly, and estimations may no longer be monotonic.

IV.3 Drawbacks and future work

A main drawback of this method is the response to noise. The standard deviations in Figure 5 are a significant fraction of the estimate itself, even though the mean is stable. However, this method will be computationally less expensive to execute over an entire PET dataset than fitting-based methods, which may outweigh its precision for low-count cases. Shibuya et al. has previously calculated that approximately $3\text{\times}{10}^{8}\text{\,}\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{s}$ will be needed to estimate oxygen concentration to a precision of $10\text{\,}\mathrm{m}\mathrm{m}\mathrm{H}\mathrm{g}$ [3]. This would correspond to three times as many counts as the extent of Figure 5. At this point, the method has essentially reached its asymptote, minimum bias, and minimum standard deviation.

Figure 10: The estimate of the $\tau_{3}=$1\text{\,}\mathrm{n}\mathrm{s}$$ case from Figure 5, however, with $\delta=0.3$ to demonstrate how the standard deviation variability can be handled for small $\tau_{3}$ values.

To control variations in the standard deviation at low counts, a small positive number $\delta$ may be introduced, as shown in Equation 9. This results in an increased bias and decreased standard deviation seen in Figure 10. Despite the increased overall bias in the asymptotic region (37.6% error), the curve remains monotonic and asymptotic, allowing for the bias to be easily corrected. For comparison, with $1\text{\times}{10}^{3}\text{\,}\mathrm{h}$istogram counts, the standard deviation decreased from $0.62\text{\,}\mathrm{n}\mathrm{s}$ for $\delta=0$ to $0.19\text{\,}\mathrm{n}\mathrm{s}$ for $\delta=0.3$, while the percentage error increased from 42.9% to 62.4% for the two cases, respectively.

Future work on this estimation method will include its testing on full three-dimensional reconstruction, testing its computational speed with respect to current estimation techniques, and further optimization of fitting parameters, such as the lower truncation point, the $e^{-st}$ term, smoothing of the histogram, and noise handling modifications.

V Conclusion

In this report we present an analytical method to estimate the orthopositronium lifetime in PALS measurements. This method uses moments of the histogram of arrival time differences, and employs an exponential weighting to mitigate numerical instability in calculation of moments from noisy data. The moment-based method was characterized in this work, and it was shown to be a stable, monotonic estimate in most cases. For cases in which the standard deviation was large, modifications to the method may be employed. This method will continue development to control noise for cases with small statistical power, and will be implemented and tested for three-dimensional PET images.

Appendix

From $f_{n}(\lambda)=\mu_{n}\{\mathrm{EXP}(t;\lambda)\}=\lambda\int_{0}^{\infty}dt\ t^{n}e^{-\lambda t}$, one can immediately derive

$$\frac{df_{n}(\lambda)}{d\lambda}=\frac{1}{\lambda}f_{n}(\lambda)-f_{n+1}(\lambda).$$ (10)

Using Equation 5 in the above equation then yields

$$f_{n+1}(\lambda)=\frac{1}{\lambda}\frac{n!}{\lambda^{n}}+n\frac{n!}{\lambda^{n+1}}=\frac{(n+1)!}{\lambda^{n+1}}.$$ (11)

When $n=0$, Equation 5 also correctly yields $f_{0}(\lambda)=1$. Therefore, by induction Equation 5 is valid for $n\geq 0$.

The function $g(t;\lambda,\sigma)=\mathrm{EXP}(t;\lambda)\ast\mathcal{N}(t;\sigma)$ is known as the Exponential Modified Gaussian (EMG) and can be shown to equal

$$g(t;\lambda,\sigma)=\lambda\ e^{-\lambda t}\ e^{\sigma^{2}\lambda^{2}/2}\ h(t;\lambda,\sigma),$$ (12)

where

$$h(t;\lambda,\sigma)=\frac{1}{2}\left(1+\mathrm{erf}\left(\frac{t-\sigma^{2}\lambda}{\sqrt{2}\sigma}\right)\right)$$ (13)

and $\mathrm{erf}(t)$ is the error function [10]. At large $t$, $h(t;\lambda,\sigma)\approx 1$ because $\mathrm{erf}(t)\approx 1$. Therefore,

$$e^{-st}g(t;\lambda,\sigma)\approx\frac{\lambda}{s+\lambda}\ e^{\sigma^{2}\lambda^{2}/2}\ \mathrm{EXP}(t;s+\lambda).$$ (14)

Note that, in evaluating $\mu_{n}\{f(t)\}=\int dt\ t^{n}f(t)$ the term $t^{n}$ in the integral progressively diminishes the contribution of $f(t)$ at small $t$ as $n$ increases. Therefore, for sufficiently large $n$ one can use the approximation given by Equation 14 to evaluate $\mu_{n}\{e^{-st}g(t;\lambda,\sigma)\}$, which then yields Equation 6.