Enhanced positronium lifetime imaging through two-component reconstruction in time-of-flight positron emission tomography

Positronium lifetime imaging (PLI) is a newly demonstrated technique possible with time-of-flight (TOF) positron emission tomography (PET), capable of producing an image reflecting the lifetime of the positron, more precisely ortho-positronium (o-Ps), before annihilation, in addition to the traditional uptake image of the PET tracer. Due to the limited time resolution of TOF-PET systems and the added complexities in physics and statistics, lifetime image reconstruction presents a challenge. Recently, we described a maximum-likelihood approach for PLI by considering only o-Ps. In real-world scenarios, other populations of positrons that exhibit different lifetimes also exist. This paper introduces a novel two-component model aimed at enhancing the accuracy of o-Ps lifetime images. Through simulation studies, we compare this new model with the existing single-component model and demonstrate its superior performance in accurately capturing complex lifetime distributions.

1 Introduction

Positron emission tomography (PET) is a major medical imaging modality that has been routinely used for diagnosing diseases such as cancer and Alzheimer’s disease [1]. In PET, one administrates molecules labeled with a positron-emitting isotope that targets certain signatures associated with these diseases so that the presence and/or the grade of disease can be manifested by the tissue uptake of the molecules at an appropriate time post-administration. PET imaging is based on coincidence detection of gamma-ray photons that are created when a positron released by the isotope annihilates with an electron. The released positrons do not annihilate instantly, and until recently, little attention has been paid to the processes leading to annihilation. In tissues, when sufficiently slowing down after traveling a certain distance (known as the positron range that is characteristic of the isotope), a positron can annihilate directly with an electron. As an alternative to such direct annihilation (DA), some positrons may form a short-lived meta-stable particle called positronium. There are two types of positronium: para-positronium (p-Ps) and ortho-positronium (o-Ps). The duration for which a positron exists before annihilation τ can be described probabilistically by an exponential distribution, p(τ) = λ e^−λτ, τ ≥ 0 so that the mean duration of existence, called lifetime, equals 1/λ. The three populations described above have different rate constants, with those of the DA and p-Ps populations substantially larger than that of o-Ps (or the DA and p-Ps populations have substantially shorter lifetimes than o-Ps). Significantly, the rate constant of o-Ps is sensitive to the local tissue properties surrounding o-Ps. The decay constant for DA and p-Ps also changes in the medium. However, the changes in DA and p-Ps decay constants in a medium are much less pronounced than the change in the o-Ps decay rate [2]. Positronium lifetime imaging (PLI) aims to exploit this feature, producing rate-constant image (or its inverse, the lifetime image) of o-Ps for providing diagnostic information for disease in addition to that given by the uptake of PET molecules. Readers are referred to Moskal et al. [3] for in-depth discussion of the physics underlying PLI.

As illustrated in Figure 1B, annihilation of a positron releases two gamma-ray photons with an energy of 511 keV that travel in opposite directions [4]. Based on this, PET detectors capture two 511-keV photons in coincidence to define a line of response (LOR) along which the annihilation is expected to take place. PET imaging describes the distribution of PET molecules in the body, and the resulting image is commonly referred to as an activity image [5]. In conventional PET, the detection time of the gamma rays is used only for determining the coincidence of two photons and the LOR. It does not determine which location on the LOR is the origin of the two photons, and hence of the positron if the positron range is negligible [6]. In time-of-flight (TOF) PET, the accuracy in time measurement is improved such that a certain degree of localization of the positron on the LOR is based on the difference in the detection time of two photons in coincidence. In the standard practice, an LOR is divided into multiple segments, which are often called a line of segment (LOS), and events detected at the LOR are further assigned to these segments according to their TOF values. The imaging model now relates the radioactivity distribution to the number of events detected at all the LOS of the system. TOF-PET is beneficial because it can increase the signal-to-noise ratio (SNR) of the resulting activity image. Better TOF resolution, typically reported using the coincidence resolving time (CRT) in picosecond (ps), leads to a better image SNR. Improving the CRT has been a hot topic for research. In the past decades, the CRT of the clinical TOF-PET system has been improved from approximately 600 ps to slightly better than 200 ps [7].

Figure 1

Figure 1. Two-dimensional illustration of the emission and detection of prompt γ (blue) and the annihilation pair (red). (A): emission of the prompt γ ray. (B): emission of the annihilation pair (two photons with an energy of 511 keV) each in opposite directions.

Recently, PLI was introduced [8–10] and ex vivo and in vivo positronium lifetime images were demonstrated using a novel multi-photon J-PET system [11, 12]. Contrary to the traditional PET/TOF-PET, PLI employs a non-pure positron-emitter, such as Sc-44, that emits a positron and a prompt gamma ray essentially at the same time and detects three-photon coincidences that include the prompt gamma and the two 511-keV annihilation photons, in addition to standard two-photon coincidences [13, 14]. In three-photon events, the detection time of the prompt gamma (depicted in blue in Figure 1A) can provide the “start” reference for measuring the elapse time of the positron before annihilation (depicted in red in Figure 1B) [2, 15]. As described above, the lifetime distribution of positrons contains three populations, and the interest in PLI is to obtain the lifetime image of the o-Ps population. In this paper, we present a two-component model that includes a slow-decay component representing the o-Ps population and a fast-decay component that represents the DA and p-Ps populations combined. The lifetime of o-Ps is very sensitive to the size of the voids (free volume between the atoms) and the concentration of molecules such as oxygen molecules [16]. There are pieces of evidence that o-Ps may provide information about disease progression in an initial stage [17, 18]. It has been shown in recent studies that the o-Ps lifetime in healthy adipose tissue differs from the o-Ps lifetime in cardiac myxoma tumors [19].

The FWHM in centimeter (cm) of the localization uncertainty illustrated in Figure 1 equals 0.015 × CRT when the CRT is expressed in ps. The best CRT of the state-of-the-art TOF-PET systems is approximately 200 ps [20, 21], corresponding to a localization uncertainty of 3 cm. Therefore, an LOS sees a mixture of events in a substantial region within which the o-Ps lifetime generally varies. Prior research studies on PLI by Qi and Huang [22]; Huang et al. [23]; Chen et al. [24]; and Shopa and Dulski [25, 26] have presented probabilistic formulations that describe the PLI data acquired by TOF-PET systems having a finite CRT, allowing accurate reconstruction of the rate-constant image. These works predominantly focused on the lifetime of o-Ps and exclusively modeled only o-Ps decays in the simulated data. In this paper, we consider more realistic simulation data by including all three positron populations.

2 Methods

2.1 Data simulation

PLI events of the o-Ps, DA, and p-Ps populations were simulated according to their respective population weights and rate-constant images by using the Monte Carlo methods and scanner configurations previously described in Huang et al. [23]. Briefly, we considered two-dimensional (2D) imaging and only triple-coincidence events. The scanner contained 364 identical detectors on a diameter of 57.2 cm. The activity image was first scaled so that its summed intensity equaled the desired number of events to produce. The number of decay in an image pixel was then drawn from a Poisson distribution whose mean equaled the scaled pixel intensity. For each decay, the prompt gamma and annihilation photons were simulated to emit in two random angles independently drawn from a uniform distribution over [0, 2π). The detectors seeing these radiations were then determined accordingly. As discussed above, the elapse time between the prompt gamma and annihilation photons was sampled from a three-component exponential mixture model using the population weights and rate constants specified for the pixel in which the decay was located. To account for uncertainty in time measurement, every detection time was perturbed by a random number drawn from a Gaussian distribution having zero mean and standard deviation (SD) σ₁. For a specified coincidence resolving time (CRT) for the scanner, the TOF bin size was chosen to be half of the CRT and ${\sigma }_1={\sigma }_t/\sqrt{2}$, where σ_t = CRT/2.35 was the SD value corresponding to the CRT. Attenuation, scattering, random events, and position range were not simulated.

PLI events of the three positron populations were generated using representative population weights (DA: p-Ps: o-Ps ≈ 0.6: 0.1: 0.3) and rate constants (λ_DA ≈ 2.5 ns⁻¹; λ_p-Ps ≈ 8.0 ns⁻¹; λ_o-Ps ≈ 0.5 ns⁻¹), as reported in the literature [10]. Figure 2 displays the ground truth of the activity image and rate-constant images and population-weight images for each component. The activity image contains two discs on an elliptical background. Pixels in the two discs have two times the activity in the elliptical background. The o-Ps rate-constant image consists of the same-sized discs and background ellipse. The two discs have different λ values (0.4 and 0.6 ns⁻¹) from the background ellipse (0.5 ns⁻¹). Outside the ellipse, since no annihilations occur, λ is not measured. Therefore, we arbitrarily set λ = 0 ns⁻¹ for this region. As mentioned above, the DA and p-Ps rate constants are not sensitive to the environment. Therefore, we consider uniform rate-constant images for them; both consisting of only the background ellipse having λ values of 2.5 ns⁻¹ (for DA) and 8.0 ns⁻¹ (for p-Ps). Similarly, the population weights are not sensitive to the environment and therefore uniform population-weight images are used, with values of 0.3, 0.6, and 0.1 for o-Ps, DA, and p-Ps, respectively. The elliptical region where λ > 0 in the rate-constant images is referred to as the region of interest (ROI). All the images are discretized into square pixels of 3.27 × 3.27 mm², with 41 × 41 pixels.

Figure 2

Figure 2. (A) Example of the ground truth of the activity image, with an average of 500 decays per pixel. (B) Ground truth of the rate-constant and population-weight images for o-Ps (OP), direct annihilation (DA), and p-Ps (PP).

Since p-Ps has a low population weight and also decays substantially faster than others, its presence in the PLI data is not significant, and we stipulate that it is sufficient to consider a PLI model that includes two components, including a fast-decay component representing the DA population and a slow-decay component representing the o-Ps population. We evaluate the performance of using this new two-component model for reconstructing the o-Ps rate-constant image.

We implemented the proposed method using the combinations of the following simulation settings.

1) Activity image used for rate-constant image reconstruction:

i) true activity image (referred to as f);

ii) estimated activity image (referred to as $\widehat{\mathit{f}}$).

2) Event size, calculated as the product of the number of pixels in the ROI (627 in the phantom used), and a multiplier of choice, which represents the average number of PLI events per pixel.

i) 627 × 100;

ii) 627 × 200;

iii) 627 × 500; and

iv) 627 × 1,000;

3) TOF resolution expressed in SD σ_t, equal to CRT/2.35

i) 0.085 ns, representing a 200-ps CRT, approximately the best TOF currently achievable

ii) 0.16 ns, representing a 377-ps TOF PET system; and

iii) 0.242 ns, representing a 570-ps TOF PET system.

For each simulation setting, we conducted 20 independent runs. The final reconstructed rate-constant images were derived by averaging the rate-constant images obtained from these 20 runs.

2.2 Two-component reconstruction model

As described by Huang et al. [23], when focusing exclusively on o-Ps decays, the elapse time measured for event k, τ_k, is modeled using an exponentially modified Gaussian (EMG) distribution, which is the convolution of an exponential distribution describing the lifetime distribution with a Gaussian distribution describing the uncertainty in measurement of elapse time [3]. In the extended model that includes both slow- and fast-decay components, τ_k follows a mixture of two EMG distributions: one representing the o-Ps population and another the DA population. As mentioned, even though the p-Ps population is included in the simulated data, in this study, we exclude the p-Ps population from the reconstruction model due to their negligible presence in the data. Our results below will demonstrate that this simplification does not noticeably compromise reconstruction accuracy.

The density of τ_k is expressed as

$w_{1,j}\cdot \mathrm{E}\mathrm{M}\mathrm{G}\left({\tau }_k;{\lambda }_{1,j},{\sigma }^2\right)+w_{2,j}\cdot \mathrm{E}\mathrm{M}\mathrm{G}\left({\tau }_k;{\lambda }_{2,j},{\sigma }^2\right),(1)$

where w_1,j + w_2,j = 1 and EMG (⋅) is the density of a EMG distribution and has the following form:

$\mathrm{E}\mathrm{M}\mathrm{G}\left(\tau ;\lambda ,{\sigma }^2\right)=\frac{1}{2}\lambda e^{-\lambda \left(x-\frac{1}{2}{\sigma }^2\lambda \right)}\left(1+\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{\tau -\lambda {\sigma }^2}{\sqrt{2}\sigma }\right)\right),(2)$

where the value for σ is determined from the CRT, which equals $(\sqrt{3}/2){\sigma }_t$. Let ${\mathcal{W}}_{N_k}=\{(c_k,{\tau }_k):k=1,\dots ,N_k\}$ denote the acquired list-mode containing N_k events where c_k and τ_k are the LOS and the measured elapse time for event k, respectively. It is shown in Huang et al. [23] that the log-likelihood function of ${\mathcal{W}}_{N_k}$ given the activity image f is

$\ell \left({\mathit{\lambda }}_{\mathbf{1}},{\mathit{w}}_{\mathbf{1}};{\mathit{\lambda }}_{\mathbf{2}},\mathit{f},{\mathcal{W}}_{N_k}\right)={\displaystyle \sum _{k=1}^{N_k}}\log \left\{{\displaystyle \sum _{j=1}^{N_j}}{H}_{c_k,j}{f}_j\left[w_{1,j}\cdot \mathrm{E}\mathrm{M}\mathrm{G}\left({\tau }_k;{\lambda }_{1,j},{\sigma }^2\right)\right.+\left.w_{2,j}\cdot \mathrm{E}\mathrm{M}\mathrm{G}\left({\tau }_k;{\lambda }_{2,j},{\sigma }^2\right)\right]\right\}.(3)$

Therefore, given an estimate for the activity image $\widehat{\mathit{f}}$ and assuming that λ₂ is known, the maximum likelihood estimation (MLE) of the o-Ps rate-constant image ${\widehat{\mathit{\lambda }}}_1$ and its weight ${\widehat{\mathit{w}}}_1$ is obtained by

$\left({\widehat{\mathit{\lambda }}}_{\mathit{1}},{\widehat{\mathit{w}}}_1\right)=\arg {\max }_{\left({\mathit{\lambda }}_{\mathit{1}},{\mathit{w}}_{\mathbf{1}}\right)}\ell \left({\mathit{\lambda }}_{\mathbf{1}},{\mathit{w}}_{\mathbf{1}};{\mathit{\lambda }}_{\mathbf{2}},\widehat{\mathit{f}},{\mathcal{W}}_{N_k}\right).(4)$

The notation used here follows the conventions established in Huang et al. [23].

• k = 1, … , N_k is the index of the PLI events.

• N_k is the total number of PLI events acquired.

• c_k is the PET-TOF LOS in which the kth PLI event is detected.

• τ_k is the measured elapse time of the kth PLI event.

• ${\mathcal{W}}_{N_k}=\{(c_k,{\tau }_k):1\le k\le N_k\}$ is the set of acquired PLI events.

• λ_1,j is the rate constant of the o-Ps population in the jth pixel. λ₁ is called o-Ps rate-constant image.

• j = 1, … , N_j is the index of the pixel.

• N_j is the total number of image pixels.

• λ_2,j is the rate constant of the DA population in the jth pixel. λ₂ is called the DA rate-constant image.

• w_1,j is the population weight of o-Ps in the jth pixel. w₁ is called the o-Ps weight image.

• w_2,j is the population weight of DA positrons in the jth pixel. w₂ is called the DA weight image.

• f_j is the number of decays in the jth pixel. f is the activity image.

• σ is Gaussian error of the EMG distribution, which is set to value $\sqrt{3/2}{\sigma }_t$, as shown in [24].

• H_c,j is an element of the system matrix H that is proportional to the probability that an isotope decay occurring inside pixel j would give rise to a two-photon coincidence event at LOS c.

w_1,j and w_2,j in Eq. (1) and Eq. (2) are the proportions of PLI events originating in pixel j that come from the o-Ps and DA populations, respectively. Evidently, w_1,j + w_2,j = 1, yielding w₂ = 1 − w₁ in Eq. (3).

In this paper, we assume that the DA rate-constant image λ₂ is known (λ₂ = 2.5 ns⁻¹) and focus on the reconstruction of the o-Ps rate-constant image. However, the proposed method can be readily extended to the scenario where the DA and o-Ps rate-constant images are both unknown. This will be investigated in future work.

2.3 Maximum likelihood estimation and performance evaluation

The MLE of f, denoted as $\widehat{\mathit{f}}$, was obtained using 50 iterations of the expectation–maximization algorithm, Lange and Carson [27], implemented in-house in Python. A total of 50 iterations were empirically determined as sufficient to reach convergence. No explicit noise regularization schemes were employed. The MLE of λ₁, denoted as $\widehat{\mathit{\lambda }}$, and its weight w was obtained by solving Eq. (4) using the limited-memory Broyden–Fletcher–Goldfarb–Shanno bound (L-BFGS-B) algorithm [28] that is available from scipy.optimize in Python, due to its ability to deal with large numbers of variables and for imposing the positivity condition for the solution (implemented using the bound constraints).

To evaluate the reconstruction performance, for each simulation, we calculated the normalized mean square error (NMSE) and the contrast recovery coefficient (CRC, NEMA Standards Publication NU 2-2001 [29]). NMSE measures the mean squared error of the estimated image in comparison with the true image, normalized by the squared Euclidean norm of the true image. CRC measures the degree of recovery of the true contrast of two circular regions compared to the background. Lower values of NMSE and higher values of CRC are indicative of better performance. Additionally, we proposed using a novel performance metric the standardized absolute log ratio (SALR) that is defined as

${\text{SALR}}_p=\frac{\left|\log \left(\frac{{\overline{\widehat{\mathit{\lambda }}}}_{1,p}}{{\overline{\widehat{\mathit{\lambda }}}}_b}\right)\right|}{\mathrm{S}\mathrm{D}\left({\widehat{\mathit{\lambda }}}_b\right)/{\overline{\widehat{\mathit{\lambda }}}}_b},p=\mathrm{1,2},(5)$

where SD stands for standard deviation, $\bar{\cdot }$ means the average, p = 1, 2 represents the left and right circle in the phantom, and the subscript b represents the background in the phantom. In Eq. (5), the denominator represents the background variability of the reconstructed image due to noise, as described in Raczyński et al. [30]. The SALR assesses the contrast of two circular regions relative to background variability. We calculate an overall SALR by averaging SALR₁ and SALR₂. Unlike other metrics, the SALR does not require knowledge of the ground-truth image. Similar to CRC, higher SALR values indicate better performance. Additional details on the computation of NMSE and CRC can be found in the Supplementary Material.

3 Results

3.1 Two-component decay model outperformed the single-component decay model

We first implemented the two-component model with the assumption that w₁ is known and used its true value. Thus, the only unknown variable was the o-Ps rate-constant image. The true w₁ = 0.3, w₂ = 0.7, and λ₂ = 2.5 ns⁻¹ were used for obtaining the MLE of λ₁. In addition, the estimated activity image $\widehat{\mathit{f}}$ was used.

Figure 3 compares o-Ps rate-constant images obtained using the two-component model versus using the previous single-component model. Figure 3A1, A2 show the averaged reconstructed image from 20 simulation runs (called mean images below). Figure 3B1, B2 plot the pixel values across the center row of the true and average reconstructed image (called horizontal profiles below). The gray shaded regions in these horizontal profiles represent the ±1 SD range of the estimated rate-constant image, calculated from the 20 simulation runs. Visually, the mean image of the two-component model agrees with the ground truth, while that of the single-component model barely resolves the two circles in the phantom. The horizontal profile of the two-component model also quantitatively agrees with the ground truth, while the profile of the single-component model is almost flat and shows substantial bias in the circles and in the background. The deviation is particularly pronounced for the right circle, which has a rate constant of 0.4 ns⁻¹. We believe that the result of the single-component model is biased toward the DA population that is not accounted for. This DA population has a large rate constant of 2.5 ns⁻¹; therefore, it generally leads to a positive bias in the reconstructed o-Ps rate-constant image. Compared to the single-component model, the two-component model yields greater SDs. However, they are small relative to the contrast of the circles, indicating stability of reconstruction.

Figure 3

Figure 3. Average of 20 reconstructed o-Ps rate-constant images and their horizontal profiles across the center of image obtained from 20 simulation runs. The gray shaded areas in the horizontal profile plots are the ±1 SD range computed from the 20 images. (A1, A2) are results obtained by the single-component model. (B1, B2) are results obtained by the two-component model. For the two-component model, the true w₁, w₂, and λ₂ were used. For both models, estimated activity image $\widehat{\mathit{f}}$ was used. Note that, as described in the text, data were simulated using three positron populations.

Figure 4 further shows box plots of the three performance metrics–NMSE, CRC, and SALR, obtained from 20 reconstructed images using the single-component and two-component models. In addition to using the true activity image f, we also used the estimated activity image $\widehat{\mathit{f}}$ obtained using the OSEM algorithm for the realistic scenario in which the true f is unavailable. The two-component model demonstrated significantly superior performance in terms of all three metrics. On the other hand, no significant difference in these performance metrics was observed between using the estimated $\widehat{\mathit{f}}$ (top row) and true f (bottom row).

Figure 4

Figure 4. Performance of reconstruction using the single-component model and the two-component model evaluated by three performance metrics. For the two-component model, the true w₁, w₂, and λ₂ were used. (A) Using the estimated $\widehat{\mathit{f}}$. (B) Using the true $\mathit{f}$.

3.2 Two-component decay model when weights are unknown

Results shown above were obtained by assuming that the o-Ps weight image w₁ is exactly known. We relaxed this assumption to allow unknown w₁ to be estimated together with λ₁ by MLE. However, we assumed a uniform weight image w₁ whose shape is known (the background ellipse), but the value is not. We believe that this is a reasonable assumption because the external boundaries of the rate-constant and weight images can be readily determined from the activity image.

We extended our investigation to evaluate the performance of the two-component model across various PLI event sizes and TOF resolutions. The previous study by Huang et al. [23] used 1 million simulated PLI events, which is equivalent to over 1,500 events per pixel. We conducted experiments with event sizes ranging from 100 to 1,000 events per pixel to evaluate the practicality and efficiency of the proposed model. Additionally, we examined the effect of TOF resolution on model performance.

Figure 5A shows the rate-constant images of a single simulation run obtained using the two-component model, allowing unknown o-Ps population weight. Various numbers of events per pixel, ranging from 100, 200, 500 to 1,000, were examined. The images in Figure 5A show stability when the number of events per pixel is decreased from 1,000 to 100. In all cases, the two circles can be readily identified.

Figure 5

Figure 5. Reconstructed rate-constant images and horizontal profiles across the center of the reconstructed images for various numbers of events per pixel.(A) Rate-constant images from a single simulation run. (B) Horizontal profiles of true and reconstructed rate-constant image.

Figure 5B shows the horizontal profile of the mean images, along with the ±1 SD range, for the four event sizes examined. The horizontal profile for 100-event-per-pixel case agrees with the truth, with a larger SD in comparison with the 1000-event-per-pixel case. Compared to Figure 3 (B2) that is obtained using true w₁, a small positive bias can be observed. These results demonstrated that the two-component method has potential to produce quantitatively an accurate o-Ps rate-constant image with a limited number of events.

Figure 6 shows the three performance metrics when the number of pixels per pixel varies. Again, the true f or estimated $\widehat{\mathit{f}}$ was used. Generally, reduction in the event size degrade the reconstruction performance and introduces greater variability for all three metrics. However, the degradation in the performance metrics appears mild, except for the 100-event-per-pixel case. Notably, unlike NMSR and CRC, the SALR demonstrates a clear monotonically increasing trend as the number of events per pixel increases. In addition, it has a very constant range of variability. Such results highlight the interpretability of SALR for evaluation of image quality and image reconstruction methods.

Figure 6

Figure 6. Box plots of NMSE, CRC, and SALR for various numbers of events per pixel, using true and estimated f.

In addition to evaluating the performance using the global CRC, SALR, and NMSE calculated from the entire reconstructed images, we also calculated the SD and bias of the three regions of interest in the reconstructed image (left circle, background, and right circle) to assess the local performance properties. The results are shown in Supplementary Figure S1 in Supplementary Material. In general, the SD decreases with the reduction of the event sizes, and bias does not decrease drastically, except for the event per pixel of 100. Details of the calculation of the bias and SD are provided in the Supplementary Material.

Furthermore, we examined three TOF resolutions, i.e., σ_t = 0.085 ns, 0.16 ns, and 0.242 ns, investigating how the reconstruction performance is affected. Figure 7A shows that the new two-component model captures the two circles even when σ_t is as large as 0.242 ns (570 ps CRT). An increasing positive bias of the reconstructed $\widehat{\mathit{\lambda }}$ can be observed in Figure 7 as TOF resolution increases, which indicates that under higher TOF resolution, the reconstructed $\widehat{\mathit{\lambda }}$ is more affected by the slow-decay component, which has a higher decay-rate constant. Figures 7A, B also show slightly increasing noise and variability in $\widehat{\mathit{\lambda }}$ with increasing TOF resolution.

Figure 7

Figure 7. Reconstructed rate-constant images and horizontal profiles across the center of the reconstructed images for various TOF resolutions. (A) Rate-constant images from a single simulation run. (B) Horizontal profiles of true and reconstructed rate-constant image.

In general, Figure 8 demonstrates a decline in the performance of three metrics, with the exception of an improvement in CRC, as TOF resolution increases from 0.16 ns to 0.242 ns. This can be explained by the fact that CRC evaluates the degree of recovery in the comparison between two circular regions and the background of the reconstructed image. A higher CRC does not necessarily indicate a higher-quality image, as it does not account for the quantitative accuracy and variability of the reconstructed image.

Figure 8

Figure 8. Box plots of CRC, NMSE, and SALR for various TOF resolutions, using true and estimated f.

3.3 Two-component decay model when weights are unknown and non-uniform

Last, we allowed non-uniform weight image w₁ but with the known shape. This assumption greatly increased the number of variables and introduced great complexity in estimation. Therefore, a larger event size of 2,000 events per pixel was used. Again, 20 simulation runs were used to compute the mean and SD images.

Figure 9A shows that the mean o-Ps rate-constant image maintained a high contrast of the two circles. The horizontal profile in Figure 9B again agreed well with the ground truth, and the ±1 SD range is reasonably small relative to the contrast of the circles.

Figure 9

Figure 9. Rate-constant image obtained by the two-component model, allowing a non-uniform o-Ps weight image. (A) Mean image from 20 simulation runs. (B) Horizontal profile across the center of the image and the ±1 SD range.

Figure 10 shows the estimated weight image ${\widehat{\mathit{w}}}_1$ and its absolute difference with the ground truth w₁. Figure 10A illustrates a high degree of accuracy of the estimated weight image as it closely approximates the true value of 0.3 across the elliptical background. The absolute difference image in Figure 10B further elucidates this precision. The results in Figures 9, 10 highlight that despite the greatly increased number of variables in estimation, the proposed two-component model with ML approach can still yield reliable solutions.

Figure 10

Figure 10. Weight reconstruction results for o-Ps using the two-component model from a single simulation run. (A) Reconstructed weight image. (B) Absolute value of the difference between the reconstructed weight image and the ground truth.

4 Discussion

This study presents a significant advancement in PLI through the development of a two-component reconstruction model in TOF-PET. By incorporating both a slow-decay and short-decay component, our model offers a more accurate description of the PLI data. The simulation studies conducted in this research clearly demonstrate the superior performance of the two-component model in capturing the intricacies of the o-Ps rate-constant image compared to traditional single-component models. This highlights the potential of our model in enhancing the accuracy of PLI, thereby contributing to more precise and informative medical imaging. The results of this study open up new possibilities for improving disease diagnosis and treatment planning, ultimately leading to better patient outcomes. Our work also lays the groundwork for further research in this area, which could lead to the development of even more advanced imaging techniques and technologies.

In this paper, the PLI data and rate-constant images are restricted to 2D cases. Further investigations are warranted to extend the proposed method to 3D settings. Additionally, the current two-component decay model treats λ₂, the slow-decay component, as known. The model can be readily extended to incorporate λ₂ as an unknown variable by providing its derivative to the optimization algorithm. However, this extension will introduce substantial computational complexities, which will be explored in future work.

Future studies should further investigate the robustness and accuracy of the model under low-dose conditions to ensure its applicability in these imaging scenarios. Additionally, extending the model to include three positron populations would allow for a more detailed characterization of the positronium decay processes for yielding further improvements in PLI. It is also of interest to add regularization to the objective function for estimating the rate-constant image.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

ZC: writing–original draft, writing–review and editing, conceptualization, formal analysis, methodology, software, and visualization. C-MK: conceptualization, methodology, software, and writing–review and editing. Hsin-Hsiung H-HH: conceptualization, funding acquisition, investigation, methodology, and writing–review and editing. LA: conceptualization, funding acquisition, investigation, methodology, and writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported in part by NSF grants DMS-1924792 and DMS-2318925 (H-HH), NIH grant R01-EB029948 (C-MK), and the United States Department of Agriculture ARZT-1361620-H22-149 (LA).

Acknowledgments

We thank Dr. Zheyuan Zhu for his help with the implementation of the proposed algorithm.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2024.1429344/full#supplementary-material

References

1. Alavi A, Werner TJ, Stȩpień EŁ, Moskal P. Unparalleled and revolutionary impact of PET imaging on research and day to day practice of medicine. Bio-Algorithms and Med-Systems (2022) 17:203–12. doi:10.1515/bams-2021-0186

CrossRef Full Text | Google Scholar

2. Bass SD, Mariazzi S, Moskal P, Stȩpień E. Colloquium: positronium physics and biomedical applications. Rev Mod Phys (2023) 95:021002. doi:10.1103/RevModPhys.95.021002

CrossRef Full Text | Google Scholar

3. Moskal P, Bozena J, Stepien E, Bass S. Positronium in medicine and biology. Nat Rev Phys (2019) 1:527–9. doi:10.1038/s42254-019-0078-7

CrossRef Full Text | Google Scholar

4. Shukla AK, Kumar U. Positron emission tomography: an overview. J Med physics/Association Med Physicists India (2006) 31:13. doi:10.4103/0971-6203.25665

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Daghighian F, Sumida R, Phelps ME. PET imaging: an overview and instrumentation. J Nucl Med Technol (1990) 18:5–13.

Google Scholar

6. Conti M. State of the art and challenges of time-of-flight PET. Physica Med (2009) 25:1–11. doi:10.1016/j.ejmp.2008.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Vandenberghe S, Mikhaylova E, D’Hoe E, Mollet P, Karp JS. Recent developments in time-of-flight PET. EJNMMI Phys (2016) 3:3–30. doi:10.1186/s40658-016-0138-3

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Moskal P. Positronium imaging. In: 2019 IEEE Nuclear Science Symposium (NSS) and Medical Imaging Conference (MIC); 26 Oct. 2 Nov. 2019; Manchester, UK (2019). p. 1–3. doi:10.1109/NSS/MIC42101.2019.9059856

CrossRef Full Text | Google Scholar

9. Moskal P, Kisielewska D, Y. Shopa R, Bura Z, Chhokar J, Curceanu C, et al. Performance assessment of the 2 γ positronium imaging with the total-body PET scanners. EJNMMI Phys (2020) 7:1–16. doi:10.1186/s40658-020-00307-w

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Moskal P, Kisielewska D, Curceanu C, Czerwiński E, Dulski K, Gajos A, et al. Feasibility study of the positronium imaging with the J-PET tomograph. Phys Med Biol (2019) 64:055017. doi:10.1088/1361-6560/aafe20

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Moskal P, Dulski K, Chug N, Curceanu C, Czerwiński E, Dadgar M, et al. Positronium imaging with the novel multiphoton PET scanner. Sci Adv (2021) 7:eabh4394. doi:10.1126/sciadv.abh4394

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Moskal P, Baran J, Bass S, Choiński J, Chug N, Curceanu C, et al. First positronium image of the human brain in vivo. medRxiv (2024). doi:10.1101/2024.02.01.23299028

CrossRef Full Text | Google Scholar

13. Matulewicz T. Radioactive nuclei for β+ γ PET and theranostics: selected candidates. Bio-Algorithms and Med-Systems (2021) 17:235–9. doi:10.1515/bams-2021-0142

CrossRef Full Text | Google Scholar

14. Choiński J, Łyczko M. Prospects for the production of radioisotopes and radiobioconjugates for theranostics. Bio-Algorithms and Med-Systems (2021) 17:241–57. doi:10.1515/bams-2021-0136

CrossRef Full Text | Google Scholar

15. Harpen M. Positronium: review of symmetry, conserved quantities and decay for the radiological physicist. Med Phys (2004) 31:57–61. doi:10.1118/1.1630494

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Moskal P, Stȩpień EŁ. Positronium as a biomarker of hypoxia. Bio-Algorithms and Med-Systems (2021) 17:311–9. doi:10.1515/bams-2021-0189

CrossRef Full Text | Google Scholar

17. Stepanov P, Selim F, Stepanov S, Bokov A, Ilyukhina O, Duplâtre G, et al. Interaction of positronium with dissolved oxygen in liquids. Phys Chem Chem Phys (2020) 22:5123–31. doi:10.1039/c9cp06105c

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Shibuya K, Saito H, Nishikido F, Takahashi M, Yamaya T. Oxygen sensing ability of positronium atom for tumor hypoxia imaging. Commun Phys (2020) 3:173. doi:10.1038/s42005-020-00440-z

CrossRef Full Text | Google Scholar

19. Moskal P, Kubicz E, Grudzień G, Czerwiński E, Dulski K, Leszczyński B, et al. Developing a novel positronium biomarker for cardiac myxoma imaging. EJNMMI Phys (2023) 10:22. doi:10.1186/s40658-023-00543-w

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Spencer BA, Berg E, Schmall JP, Omidvari N, Leung EK, Abdelhafez YG, et al. Performance evaluation of the uEXPLORER total-body PET/CT scanner based on NEMA NU 2-2018 with additional tests to characterize PET scanners with a long axial field of view. J Nucl Med : official Publ Soc Nucl Med (2021) 62:861–70. doi:10.2967/jnumed.120.250597

CrossRef Full Text | Google Scholar

21. Alberts I, Hünermund J-N, Prenosil G, Mingels C, Bohn K, Viscione M, et al. Clinical performance of long axial field of view PET/CT: a head-to-head intra-individual comparison of the biograph vision quadra with the Biograph Vision PET/CT. Eur J Nucl Med Mol Imaging (2021) 48:2395–404. doi:10.1007/s00259-021-05282-7

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Qi J, Huang B. Positronium lifetime image reconstruction for TOF PET. IEEE Trans Med Imaging (2022) 41:2848–55. doi:10.1109/tmi.2022.3174561

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Huang H-H, Zhu Z, Booppasiri S, Chen Z, Pang S, Kao C-M. A statistical reconstruction algorithm for positronium lifetime imaging using time-of-flight positron emission tomography (2023). arXiv preprint arXiv:2206.06463.

Google Scholar

24. Chen Z, An L, Kao C-M, Huang H-H. The properties of the positronium lifetime image reconstruction based on maximum likelihood estimation. Bio-Algorithms and Med-Systems (2023) 19:1–8. doi:10.5604/01.3001.0054.1807

CrossRef Full Text | Google Scholar

25. Shopa RY, Dulski K. Multi-photon time-of-flight MLEM application for the positronium imaging in J-PET. Bio-Algorithms and Med-Systems (2022) 18:135–43. doi:10.2478/bioal-2022-0082

CrossRef Full Text | Google Scholar

26. Shopa RY, Dulski K. Positronium imaging in J-PET with an iterative activity reconstruction and a multi-stage fitting algorithm. Bio-Algorithms and Med-Systems (2023) 19:54–63. doi:10.5604/01.3001.0054.1826

CrossRef Full Text | Google Scholar

27. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr (1984) 8:306–16.

PubMed Abstract | Google Scholar

28. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods (2020) 17:261–72. doi:10.1038/s41592-019-0686-2

PubMed Abstract | CrossRef Full Text | Google Scholar

29. NEMA Standards Publication NU 2-2001. Performance measurements of positron emission tomographs. Rosslyn, VA: National Electrical Manufacturers Association (2001).

Google Scholar

30. Raczyński L, Wiślicki W, Klimaszewski K, Krzemień W, Kopka P, Kowalski P, et al. 3D TOF-PET image reconstruction using total variation regularization. Physica Med (2020) 80:230–42. doi:10.1016/j.ejmp.2020.10.011

CrossRef Full Text | Google Scholar