SPLIT: Statistical Positronium Lifetime Image Reconstruction via Time-Thresholding

Positron emission tomography (PET) is a widely utilized medical imaging modality that uses positron-emitting radiotracers to visualize biochemical processes in a living body. The spatiotemporal distribution of a radiotracer is estimated by detecting the coincidence photon pairs generated through positron annihilations. In human tissue, about 40% of the positrons form positroniums prior to the annihilation. The lifetime of these positroniums is influenced by the microenvironment in the tissue and could provide valuable information for better understanding of disease progression and treatment response. Currently, there are few methods available for reconstructing high-resolution lifetime images in practical applications. This paper presents an efficient statistical image reconstruction method for positronium lifetime imaging (PLI). We also analyze the random triple-coincidence events in PLI and propose a correction method for random events, which is essential for real applications. Both simulation and experimental studies demonstrate that the proposed method can produce lifetime images with high numerical accuracy, low variance, and resolution comparable to that of the activity images generated by a PET scanner with currently available time-of-flight resolution.

and identification of molecular-level activities within a living body.This technology is widely utilized in the fields of oncology [1], neurology [2], and cardiology [3].The current approaches of PET concentrate on measuring the concentration of radioactive tracers.This is accomplished by using PET scanners to detect coincident photon pairs produced from positron annihilations.The resulting activity image is then generated by maximizing the likelihood of the detected events, which is governed by Poisson statistics.
Current PET scanners completely ignore the life history of positrons before annihilations.However, due to the presence of the positronium, a bound state of an electron and positron, the history of positrons before annihilations is also worth investigating.The positroniums are formed by a thermalized positron and an electron, with an occurrence up to 40% in human tissue [4].There are two kinds of positroniums: parapositronium (p-Ps) and ortho-positronium (o-Ps).The p-Ps has a relatively short lifetime of 0.125 ns [5] in vacuum before annihilating into two 511-keV photons; however, the lifetime of o-Ps varies greatly from 142 ns [6] in vacuum to several ns in tissues [7], [8], and the o-Ps annihilates into three photons in vacuum but predominantly two photons in tissues [9].This variability of o-Ps behavior is due to two effects-pick-off annihilation and spin-exchange interaction [10]-which react much faster than the three-photon annihilation and thus become dominant in materials.These two effects are both due to the interaction between a o-Ps and its surrounding microenvironment: pick-off annihilation occurs when the positron of the o-Ps annihilates with a foreign electron; and spin-exchange is induced when the surrounding molecules possess unpaired electrons.Therefore, the o-Ps lifetime is dependent on the size of intermolecular voids and the concentration of bio-active molecules in biological materials [11].Previous studies [12], [13] used positron annihilation lifetime spectroscopy (PALS) to measure the o-Ps lifetime in water samples with different O 2 concentrations and found that the inverse of lifetime increases linearly with increased dissolved oxygen concentration (pO 2 ).This result indicates that o-Ps lifetime would be an indicator of hypoxia and the ability of identifying hypoxia regions noninvasively in vivo would enable more effective treatments [14], [15], [16].The first positronium lifetime images of spatially separated cardiac myxoma samples and adipose samples showed that the o-Ps lifetimes in cardiac myxoma tissue were 1.950 ± 0.019 ns and 1.874 ± 0.020 ns for two subjects, respectively, and the o-Ps lifetimes in adipose were 2.645 ± 0.027 ns and 2.581 ± 0.030 ns for the same two subjects, respectively [8].
To measure the lifetime, positron emitters with a prompt gamma emission, which provides a start signal, are needed, such as 44 Sc, 22 Na, and 124 I [17].Instead of the conventional 511-keV coincidences in PET, triple-coincidences (tri-coincidences) each consisting of one prompt gamma and two 511-keV photons are detected, with the time difference between the detection of 511-keV photon pair and prompt gamma being a random measurement of lifetime.Due to the lower efficiency of detecting tri-coincidences as compared to the 511-keV pairs, positronium lifetime imaging poses the first challenge on detector sensitivity [18].Fortunately, this problem has been addressed by the advent of long axial field-of-view (FOV) PET scanners, such as the EXPLORER scanner [19], PennPET Explorer [20], Siemens Vision Quadra [21], and J-PET scanner under construction [22].Computer simulation shows the sensitivity of tri-coincidence detection on the EXPLORER scanner is even higher than the detection of 511-keV pairs on current whole-body PET scanners [18], [23].
This paper aims to address the second challenge: the lack of efficient image reconstruction methods for positronium lifetime imaging.We focus on imaging the lifetime using the 2-photon annihilation of o-Ps, rather than the 3-photon annihilation [9], which is extremely rare in biological tissue and methods for localizing the annihilation point using 3 photons already exist [24], [25].We also limit our analysis to lifetime imaging using radionuclides emitting prompt gammas at a higher energy than 511 keV that can be separated from annihilation photons by energy thresholding.Most existing studies based on the 2-photon annihilation focus on imaging well-separated small samples [8], [23], [26], where one can use time-of-flight (TOF) information to localize and thus differentiate different sources spatially, and the lifetime is then estimated by fitting the histogram of the lifetime measurements using a lifetime model.As a result, the spatial resolution is limited by the TOF resolution.The TOF resolutions (full-width-at-half-maximum (FWHM) of coincidence resolving time (CRT)) of modern PET scanners are 200-500 ps, corresponding to 30-75 mm uncertainties in images and are far from sufficient for imaging large heterogeneous objects.Furthermore, due to the poor spatial resolution, events from regions of different lifetimes are likely to be mixed, which would make the lifetime model more complicated by adding more o-Ps lifetime components.To obtain high-resolution lifetime images using PET scanners with a currently available TOF resolution, we have previously derived a penalized maximum likelihood (PML) method for positronium lifetime image reconstruction [27].As a statistical reconstruction method, it can recover high-resolution and low-variance images.However, it has several limitations: (1) its high computational cost prevents its usage in fully 3D reconstruction; (2) it assumes a single exponential distribution, which is too simple to characterize the positron annihilation in real situations; (3) it does not have a random events correction method, which is essential for real applications of PLI.In this paper, we propose a new lifetime image reconstruction method that is more computationally efficient than the previous PML method.We call it statistical positronium lifetime image reconstruction via time-thresholding (SPLIT).This method features a straightforward implementation leveraging the existing 3D PET image reconstruction algorithms and a new random events correction method.

A. PLI Event Model
In this paper, a PLI event is defined as a tri-coincidence consisting of two annihilation photons and a prompt gamma.It is mathematically represented by a TOF sinogram bin i k determined by the line of response (LOR) and the TOF bin of the coincidence event and a lifetime measurement τ k which is the estimated time difference between the production of a prompt gamma and the subsequent annihilation after travel time correction (see Section II-C).We first focus on true PLI events consisting of three photons resulted from the same decay.The analysis for random PLI events will be presented later.For true PLI events, the lifetime measurement τ follows a distribution modeled as a summation of multiple exponential functions convolved with a Gaussian function [28], [29]: where g (τ ) is a gaussian function modelling the detector timing response and u (t) is the unit step function.The subscript l ∈ {o, p, d} denotes an annihilation pathway corresponding to o-Ps, p-Ps and direct annihilation, respectively.A l denotes the intensity of the l th pathway and l∈Ê{o, p,d} A l = 1.λ l is the annihilation rate (i.e., the inverse of the lifetime).

B. SPLIT Algorithm
The core of the SPLIT method is to decouple the image reconstruction and lifetime estimation by obtaining a lifetime-encoded activity image z j (T c ) using the PLI events with lifetime measurements in the window [T 1 , T c ] , where T 1 is a constant lower threshold and T c is a variable upper threshold.For an event originated in voxel j, the probability for it to be retained after the thresholding is where p τ | λ j , A j is the lifetime distribution defined in (1).Usually, T 1 is set to be low enough to include all PLI events with lifetime measurements less than T c , but it can also be set to zero or even a higher value.Let x j = q j x 0 j , where x 0 j is the total activity in voxel j and q j is the sensitivity of detecting a prompt gamma from voxel j by any detector in the prompt gamma energy window.We introduce a lifetime-encoded image z j (T c ) = P T c ; λ j , A j x j and obtain the following relationship between z j (T c ) and the expected number of the retained PLI events under threshold T c : where ḡi (T c ) is the expectation of the number of retained events in TOF sinogram bin i within the window [T 1 , T c ]; H i j is the (i, j) th element of the standard TOF PET system matrix, denoting the probability of detecting a 511-keV photon pair in TOF sinogram bin i from a positron annihilation in voxel j.It is worth noting that the above model in (3) uses the standard PET system matrix and does not require the calculation of the prompt gamma detection efficiency q j .The advantage of using the standard PET system matrix is that the activity image z (T c ) can be reconstructed by any existing PET image reconstruction algorithm, such as the ordered-subset expectation maximization (OSEM) algorithm.
To estimate the lifetime image, a series of lifetime-encoded images are reconstructed using a set of thresholds The o-Ps lifetime in voxel j is obtained by minimizing the squared error shown below In practice, it is neither feasible nor necessary to fit the three lifetimes and two effective intensities for each voxel because the lifetimes of p-Ps and direct annihilation are very short (<0.5 ns) and almost not affected by the surrounding materials [30].Thus, we first analyze the lifetime spectrum of all collected PLI events to estimate these two short lifetimes prior to the SPLIT reconstruction and then fix their values in the subsequent threshold-activity fitting to estimate o-Ps lifetime λ o, j , and two effective intensities A o, j , A p, j within each voxel.In the rest of this paper, we will use θ j = [λ o, j , A o, j , A p, j ] to denote the parameters of interest in the SPLIT reconstruction.

C. Tri-Coincidence Travel Time Correction
Since the photons from the same decay are likely to travel different distances before being detected, it is necessary to estimate their emission times for an accurate lifetime measurement.This is a pre-reconstruction correction for lifetime measurements only.First, an annihilation point is inferred as the most likely position along the LOR using the TOF information of the two 511-keV photons.Then the emission time of each photon is estimated by subtracting the travel time between the estimated annihilation point and its detection point inside the detector ring from its detection time.Finally, the lifetime measurement of each PLI event is computed as τ = t 511 − t P G , where t 511 and t P G are the estimated times of the positron annihilation and the prompt gamma emission, respectively.Note that this correction is only performed for each lifetime measurement, but the TOF information used in the activity image reconstruction is unaltered.
Because of the uncertainty in the annihilation point estimation, the travel time correction is not perfect and results in additional time blur to the lifetime measurements.This effect can be lumped into the detector timing response g (τ ) in ( 1).As a result, the FWHM of g (τ ) is different from the TOF resolution measured at 511 keV.In this paper, we estimate the FWHM of g (τ ) in the pre-reconstruction fitting together with the lifetimes of direct and p-Ps annihilations.

D. Random Events Estimation
We developed a software-based method to group the tri-coincidences from a data stream consisting of singles data.We focus our attention on radionuclides emitting prompt gammas at a higher energy than 511 keV and the prompt gammas can be separated from annihilation photons by energy thresholding.The grouping algorithm traverses through every 511-keV single event and uses it as a reference event, from which the other 511-keV photon and the prompt gamma are searched within their respective time windows as shown in Fig. 1.These two time windows are collectively referred to as the prompt window if neither of them is delayed.We used take-all-goods policy to retrieve every eligible tri-coincidence in the prompt window.
The random events in the prompt window are more complex than those in standard PET because there are three types of random tri-coincidences in the prompt time window depending on the relationship among the three single events as shown in Table I.Type I randoms are formed by a pair of 511-keV photons from the same annihilation with an unrelated prompt gamma; Type II randoms are formed by a 511-keV photon and a prompt gamma from the same decay with an unrelated 511-keV photon; Type III randoms are formed by photons from three different decays.
To estimate the number of random tri-coincidences, we use four configurations of delayed time windows, with each designed to contain one or two types of random Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I SUMMARY OF POSSIBLE RELATIONSHIPS BETWEEN EACH TWO SINGLE EVENTS OF A TRI-COINCIDENCE. 'T' MEANS THE TWO EVENTS ARE FROM THE SAME DECAY AND 'R' MEANS THEY ARE FROM TWO DIFFERENT DECAYS (I.E., RANDOM). THE 1 ST 511-KEV SINGLE EVENT SERVES AS A REFERENCE SINGLE AND IS ALWAYS PRIOR
TO THE 2 ND 511-KEV SINGLE EVENT tri-coincidences (Fig. 1).For tri-coincidences yielded in the first two delayed windows, the two 511-keV single events are unrelated to each other due to the large delay; in the third delayed window, all three singles are unrelated to each other; in the fourth delayed window, the prompt gamma is unrelated to the 511-keV singles.The detected types of random events of each delayed window are listed on the right in Fig. 1.Since type III randoms are present in every delayed window, the number of events in delayed window #3 needs to be subtracted from the number of events in the other delayed windows to estimate type I, II.a, and II.b randoms.Similar to conventional PET coincidence data, we refer to the tri-coincidences formed by the prompt time window as prompt events and the tri-coincidences formed by the delayed time windows as delayed events.When calculating the lifetime measurements for delayed events, we first take out the time offset (50 ns or 150 ns) and then compute the lifetime measurement of each delayed event in the same way as a prompt event, including travel time correction.

E. Implementation of SPLIT Method With Random Events Correction
First, we noticed that type II and III random events are similar to random coincidences in conventional PET in that they consist of two unrelated 511-keV photons and cannot be reconstructed into a valid activity image.In comparison, the two 511-keV photons in a type I random event are from the same positron annihilation and they can be mapped to a valid activity image.Considering these properties of random events, we propose to correct for type II and III randoms in the data domain during the image reconstruction and then to model type I randoms in the fitting of the threshold-activity curves after the OSEM reconstruction.
With the consideration of randoms, the forward projection model in (3) becomes where the new lifetime-encoded image z T m c contains both the true activity and type I random events and r m i models types II and III random events in TOF sinogram bin i under threshold T m c .Because the two 511-keV photons in types II and III randoms are unrelated, such random events are uniformly distributed among TOF bins.Thus, r m i can be estimated by where [delayed window#k] m i denotes the number of tri-coincidence events with two 511-keV photons detected by the two detectors forming the i th TOF sinogram bin in the k th delayed window and within the lifetime window [T 1 , T m c ].The factor of 2 in the denominator arises from the fact that the effective coincidence window is from −T 511 to +T 511 .
Because the prompt gammas in type I randoms are not related to the 511-keV photon pairs, the lifetime measurement of the type I random events follows a uniform distribution.Therefore, the new lifetime-encoded image z (T c ) can be expressed as where B j is the rate of type I random events with the 511-keV photon pairs originated from voxel j.When T c is far greater than the o-Ps lifetime in tissue, we have P T c ; θ j ≈ 1 and z j (T c ) becomes a linear function of T c , i.e., This linear range provides a simple way to estimate x j and B j using linear regression.In this paper, we first fit a linear function to z j (T c ) for T c between 50 and 100 ns and the resulting slope is B j and the intercept ( x j − B j T 1 ).Note that the estimates x j and B j are independent of the lifetime parameters.Then we used thresholds lower than 20 ns (nonlinear range) to fit θ j by minimizing the mean squared error between the reconstructed lifetime-encoded activity curve and their expected values voxel-by-voxel: where ẑ j (T m c ) denotes the reconstructed lifetime-encoded images for the thresholds and M is the number of time thresholds in the nonlinear range.The complete workflow of the SPLIT method with randoms correction is shown in Fig. 2.

F. Methods to Compare
The SPLIT method was compared with the previously published list-mode maximum likelihood (ML) lifetime image reconstruction method [27] and the direct TOF backprojection (direct TOF-BP) method.
The ML lifetime reconstruction method assumes a mono-exponential lifetime distribution and uses an iterative algorithm to find the maximizer of the log likelihood function of the list-mode PLI events.The direct TOF-BP method first localizes a PLI event based on the most likely position provided by the TOF information along its LOR and then estimates the o-Ps lifetime in each voxel by fitting the lifetime spectrum of the events assigned to each voxel.We found that fitting the spectrum based on the weighted least squares criterion as commonly done in PALS studies resulted in poor performance because the number of PLI events per voxel is very low in this study.Here we adopted a maximum likelihood-based method to estimate the o-Ps lifetime using the following lifetime model where D j is a background intensity modeling all random events assigned to voxel j and is different from B j in (8) as the latter only models type I randoms.This model ignores the nonuniformity in the time distribution of type II random events (see Fig. 4), but the approximation is reasonable for the activity level used in this paper.The o-Ps lifetime is obtained by maximizing the log likelihood of the PLI events within each voxel where k is the index for PLI events and K is the number of events assigned to the voxel.The optimization was performed using the mle function in MATLAB.The o-Ps lifetime, the intensities of o-Ps and p-Ps, and D were fitted with the two short lifetimes being fixed.

A. Comparison With List-Mode ML Reconstruction
The SPLIT method was first compared with the list-mode ML lifetime reconstruction method using the same simulation setup described in [27].Briefly the simulation modeled a single-ring PET scanner.The phantom consists of an elliptical background with two inserts of different lifetimes, represented by an image array of 21 × 21 pixels (Fig. 3a).The phantom has a uniform activity inside the elliptical region (phantom  one in [27]).Because the list-mode ML method can only deal with lifetime spectra with one component without randoms, the simulation only included one lifetime component and the curve fitting in the SPLIT method was reduced to this simplified model as well.No random event was simulated.The ML reconstruction was performed with 100 iterations.For the SPLIT method, twelve thresholds T m c were set between 0 ns and 20 ns: (1.0, 1.2, 1.5, 1.9, 2.4, 3.0, 3.8, 5.0, 7.0, 9.0, 14.0, 20.0) ns; 20 iterations and 3 subsets were used for the OSEM reconstruction.Fifty independent identically distributed datasets were reconstructed by both the ML and the SPLIT methods to calculate the mean and standard deviation of the reconstructed lifetime images.

B. GATE Simulation
To validate the SPLIT under a more realistic 3D imaging situation, we performed a Monte Carlo simulation using GATE v8.2 [31].The PET scanner was the NeuroEXPLORER [32] with a ring diameter of 52 cm and an axial length of 49 cm.The simulated TOF resolution was 250 ps, which corresponds to a FWHM of 177 ps for the Gaussian time blur added to each single event.The energy resolution was set to 13% at 511 keV.The phantom was a 3D rodent image with an axial length of 13.6 cm.A 10-mm diameter spherical lesion Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
was inserted into the lower right flank of the body. 44Sc was chosen as the radioisotope (prompt gamma energy 1157 keV).The simulated activity distribution mimicked prostate-specific membrane antigen (PSMA) distribution with an activity ratio of 10 : 15 : 2 : 1 between the lesion, kidneys, liver, and body background as shown in Fig. 5 (a).The total activity was 2.78 MBq with an average concentration of 20 kBq/cc and the total scan time was 30 min.
The lifetime was simulated after the GATE simulation by retrospectively modifying the time tags of the prompt-gamma events according to the origin of each event and the corresponding o-Ps lifetime at the original location.The o-Ps lifetime was set to 2.0 ns inside the lesion and 2.5 ns elsewhere inside the body (Fig. 5 (c)).The p-Ps and direct annihilation lifetimes were 0.125 and 0.4 ns, respectively.The intensities of o-Ps, p-Ps, and direct annihilation were set to 30%, 10%, and 60%.In this way, a list of singles data of the PLI scan was ready for the grouping of tri-coincidences.The energy window for 511-keV photons was [430 keV, 650 keV] and for the prompt gammas, the lower energy threshold was 700 keV and there was no higher energy threshold.The width of the 511-keV time window was T 511 = 1 ns.The prompt gamma window was [−100, 10] ns.Twenty-three values of the threshold T m c were set in the range of 1-20 ns: (1.00, 1.10, 1.20, 1.35, 1.50, 1.70, 1.90, 2.15, 2.40, 2.70, 3.00, 3.40, 3.80, 4.40, 5.00, 6.00, 7.00, 8.00, 9.00, 11.00, 14.00, 17.00, 20.00) ns and 6 were set in 50-100 ns (50, 60, 70, 80, 90, 100) ns; the lower bound T 1 was set to −1 ns.The tri-coincidence data were stored in list mode with a TOF bin width of 39 ps.
The lifetime-encoded images were reconstructed by the list-mode TOF OSEM using 3 subsets and 2 iterations.Randoms were corrected using the method described in Section II-E, but no scatter correction was considered in the reconstruction.The reconstruction FOV was 50 mm × 50 mm × 180 mm and the reconstruction voxel size was 0.8 mm × 0.8 mm × 1.6 mm, which is half of the detector crystal pitch.The reconstructed lifetime-encoded images were fitted using the MATLAB function fmincon for each voxel to estimate θ j = [λ o, j , A o, j , A p, j ].The lifetimes of p-Ps and direct annihilation were fixed to pre-determined values obtained by fitting the lifetime spectrum of all PLI events by PALSfit3 [29], a program used widely in PALS experiments.In this pre-reconstruction fitting, we use three lifetime components to model direct, p-Ps, and o-Ps annihilations because more than 99% of the PLI events have an o-Ps lifetime of 2.5 ns.The fitted lifetimes of p-Ps and direct annihilation were used in the subsequent fitting for the threshold-activity curve.The lifetime image was also reconstructed using the direct TOF-BP for comparison.

IV. REAL EXPERIMENTAL STUDY
The proposed SPLIT method was also applied to the real data from an ex vivo scan on the J-PET tomograph [33], [34].The tangential, radial, and axial NEMA-based spatial resolutions of the scanner were 0.64 cm ± 0.09 cm, 0.28 cm ± 0.09 cm, and 3.05 cm ± 0.03 cm, respectively, at a radial position of 1 cm from the scanner center.The TOF resolution was estimated to be 460 ps (FWHM) [8].The sensitivity for the tri-coincidence event detection in the experiment was estimated to be 0.016 cps/kBq.In this experiment, two cardiac myxoma tissue samples and two adipose tissue samples were scanned [35].These samples were extracted from two patients, with each of them providing one cardiac myxoma and one adipose sample.A 22 Na point source covered by a layer of thin Kapton foil was inserted into each tissue sample so that positrons emitted by 22 Na were able to enter the tissue.In the transaxial plane, these four samples were placed on the four corners of a square with a side length of 16.2 cm.More details of experimental setup and scanner have been published previously [8].
Since this experiment was originally designed to be reconstructed by the direct TOF-BP method, the sources were placed at a large distance with almost no interference between them.To demonstrate the spatial resolution of the SPLIT method, we also adopted a shrinking technique to virtually bring the four point sources closer to a distance of 5.4 mm while preserving the spatial and TOF resolutions of the system.Details of the shrinkage operation are described in Appendix B. No change was applied to the lifetime measurements.
Both the original data and shrunk data were reconstructed using the direct TOF-BP and SPLIT methods.Types II and III randoms were omitted because they were insignificant as compared to type I randoms at this activity level, which can be inferred from the relative intensity in the simulation study plotted in Fig. 4. The same lifetime thresholds as in the GATE simulation study were used.The reconstructed voxel size was 2.5 mm × 2.5 mm × 5 mm.Three subsets and 5 iterations were used for the OSEM reconstruction.Based on a separate PALS measurement, four lifetime components existed in this case: the intensities of o-Ps, p-Ps, and direct annihilation in the tissue sum to 90% with the annihilation in the Kapton foil making up the rest [8].The lifetimes of p-Ps and direct annihilation were fixed to 0.125 ns [30] and 0.388 ns [8], respectively, and the lifetime and intensity of the annihilation in the Kapton foil were fixed to 0.374 ns and 10% [36], which were consistent with the values used in the previous study [37].Same as the simulated rodent scan, θ j = [λ o, j , A o, j , A p, j ] were estimated in the SPLIT fitting.The reconstructed lifetime in each region was quantified for comparison.

A. Comparison With List-Mode ML Reconstruction
The lifetime images reconstructed by the list-mode ML and SPLIT are averaged across 50 realizations and the mean lifetime images are shown in Fig 3 (b,c), with the standard deviation (s.d.) images shown in Fig 3 (d,e).The mean lifetime values of the three regions of interest (ROIs) across 50 realizations and the s.d. of the ROI means are listed in Table II.Overall, the performances of the SPLIT method and list-mode ML are comparable, although the SPLIT method may yield a slightly larger bias.We also notice that the list-mode ML method results in higher pixel-wise s.d. in the right insert because it does not model the Gaussian blurring due to the TOF resolution in the likelihood function and the model mismatch is more significant in the right insert which has a shorter lifetime.In addition to the ability to utilize an  accurate lifetime model, another main advantage of SPLIT method is its computational efficiency.In this experiment with a pre-computed system matrix, the SPLIT took 11 s for the OSEM reconstruction with twelve time thresholds and 1.1 s for the curve fitting.The ML lifetime reconstruction took 940 s for 100 iterations with CPU parallelization.Both methods were run in MATLAB on an Intel Core i9-10920X 3.5 GHz CPU.The reason for this difference is mainly because the ML lifetime update needs to re-compute an image to be forward projected for every event.

B. GATE Simulation
In Fig. 4, we plot the histograms for the lifetime measurements of the events in the prompt window ("prompts"), and three types of randoms estimated by delayed windows for the simulated rodent phantom scan.We can see that types I and III random events exhibit a uniform lifetime spectrum because the prompt gamma is unrelated to the annihilation photons, whereas type II randoms have a non-uniform lifetime spectrum because the prompt gamma is related to one annihilation photon.We also plotted the lifetime histograms of the true events (identified by the eventID in GATE) and Prompts − Randoms.Clearly, the Prompts − Randoms histogram matches with the ground truth very well, indicating the proposed delayed windows estimate the random events very accurately.The total numbers of trues, prompts, and estimated type I-III random events are 89585971, 95415822, 5731026, 91457, and 3245 respectively.The relative difference between the number of trues and the number of Prompts− Randoms is only 0.005%.It is also important to notice that at this activity level (2.78 MBq total activity), type I random events are predominant.
The p-Ps and direct annihilation lifetimes determined in the pre-reconstruction fitting were 0.1275 ns and 0.4014 ns, respectively.These fitted values are very close to their ground truths: 0.125 ns and 0.4 ns for p-Ps and direct annihilation, respectively.The estimated FWHM of g (τ ) from the pre-reconstruction fitting was 238 ps, which contains uncertainties from both TOF measurement and imperfect travel time correction.The reconstructed activity and lifetime images are shown in Fig. 5.Note that the reconstructed activity image is weighted by the prompt gamma detection efficiency factor q j with the average of q j normalized to 1. Positron range was not corrected in all activity reconstructions.The average lifetimes in different ROIs are listed in Table III.Note that in both direct TOF-BP and SPLIT lifetime images, a voxel is set to zero if its activity is below 50% of the activity in the body background.The direct TOF-BP method works well in recovering the average lifetime in large regions, but there is a significant bias in the lesion.This is because the 250-ps TOF resolution translates to a spatial resolution of 37.5 mm, resulting in many events from the normal tissue being misplaced in the lesion and vice versa.Since lesion is a relatively small region with fewer events than the normal tissue, its accuracy is more susceptible to these mispositioned events.The direct TOF-BP also suffers high variance as the average counts per voxel (0.8 mm × 0.8 mm × 1.6 mm) is only around 400.
Compared to the direct TOF-BP, the proposed SPLIT method with the pre-fitted short lifetimes shows accurate mean lifetimes with well-controlled standard deviations in both normal tissue and lesion.For comparison, we also performed SPLIT fitting by fixing the two short lifetimes to their ground truth values and the resulting o-Ps lifetimes are shown in Table III.As expected, the SPLIT with the true short lifetimes has slightly lower biases and standard deviations than the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.SPLIT with fitted short lifetimes, but the difference is relatively small, indicating the pre-reconstruction fitting method is practically feasible.

C. Real Experimental
The reconstructed lifetime images of the ex vivo tissue samples are shown in Fig. 6 and their lifetime values are listed in Table IV.To obtain a reference value of the average o-Ps lifetime for each sample, the lifetime spectrum of each source was analyzed, and the resulting lifetime values are shown in Table IV.The ROI in each quadrant was determined by the region with activity greater than 10% of the peak intensity of each source.As the sources become closer, the lifetime accuracy of direct TOF-BP method starts to degrade because the fraction of mispositioned events increases in each ROI.The direct TOF-BP method also suffers from high variance.However, the lifetime variance of the direct TOF-BP method decreases slightly as the sources become closer, because the blurring increases the number of events in each voxel (the number of events per voxel is shown in Fig. 6 (a)) and the variance of the ML estimation is predominantly determined by the total counts per voxel.The SPLIT method reconstructs accurate lifetime images at the 30-fold shrinkage, where the distance between two adjacent sources is 5.4 mm.The SPLIT method is also effective in variance reduction, yielding significantly lower variances than those of the direct TOF-BP.

VI. DISCUSSION
The general idea of the proposed SPLIT method is to encode the lifetime into the intermediate activity images and decode the lifetime out of these activity images.To do this, we used time thresholds to generate a series of list-mode datasets and reconstruct the time-encoded activity curve for each voxel.There are other ways to accomplish this goal.For example, the projection data can also be created by retaining the events with lifetime greater than a threshold or between two thresholds.The subsequent fitting step needs to be modified accordingly.The advantage of using thresholds to sort the events is to maintain the Poisson statistics of the list-mode data so that a number of established activity image reconstruction methods can be applied without any modification.
The computational burden of the SPLIT method depends on the number of lifetime-encoded images.We used a total of 20-30 thresholds in this study.In general, more thresholds lead to a better estimation of the lifetime, but the extent of improvement reduces as the number of thresholds increases.A previous study on fluorescence lifetime estimation [38] suggested 10 or more time bins for a monoexponential decay model and indicated that the number of bins for a multiexponential model (similar to the one used here) would be substantially greater.It will be worth investigating the optimal number of thresholds in the SPLIT method to reach a good balance between the computation cost and image quality.In the current SPLIT procedure, we firstly fit the complete lifetime spectrum prior to the reconstruction to obtain the lifetimes of p-Ps and direct annihilation and fix them in the subsequent fitting of the threshold-activity curves.This is because these two lifetimes are fairly short and independent of the surrounding materials.This substantially improves the quality of the lifetime estimate because the current TOF resolution is on the same level of the short lifetimes, making it difficult to recover these values along with the o-Ps lifetime in the threshold-activity curve fitting.Our simulation results showed that the pre-reconstruction fitting method is effective.
The current lifetime images are estimated using the activity images from early OSEM iterations because noise significantly deteriorates the lifetime image quality if later iterations were used.It is well known that terminating the OSEM reconstruction early would incur strong correlations between voxels [39] and a degraded spatial resolution [40].A better solution would be to incorporate prior information in the reconstruction to reduce the noise.There have been many effective regularization-based methods [41] and kernel-based methods [42] to improve the quality of activity images.It will be a future task to incorporate them into the SPLIT method.
The number of random events in PLI can be a major issue in real human applications.In Appendix A, we prove that the randoms-to-trues ratio is dependent on the widths of the tri-coincidence time windows and the total activity, but only weakly dependent on the detection efficiencies of the scanner, Furthermore, the types I and II randoms-to-trues ratios increase linearly with respect to the activity level, while the type III randoms-to-trues ratio increases quadratically with respect to the activity level.To validate this theoretical prediction, we increased the total activity of the simulated rodent scan from 2.78 MBq to 55. 6 MBq and grouped all the eligible tri-coincidences in the time and energy windows without considering the support of the object (the dimension of the object was used to reject non-colinear 511-keV pairs in the results section).The randoms-to-trues ratios for type I, II, and III randoms were 7.3%, 1.7%, and 0.055%, respectively at 2.78 MBq, and 146%, 34%, and 22%, respectively at 55. 6 MBq.The increases of randoms-to-trues ratios were in good agreement with the theoretical prediction.It shows that at the normal activity level of a human scan, the number of randoms can exceed the number of trues.Therefore, the randoms correction method presented in this paper is very important in performing positronium lifetime imaging in humans.
In our simuation study, type II and III randoms were estimated using delayed windows according to (7).We did not employ any variance reduction technique because the TOF bin width (39 ps) is much smaller than T 511 (= 1ns), resulting in a much smaller variance in the randoms estimate than in the prompt data.When needed, further variance reduction is possible by increasing the width of delayed 511-keV window and other variance reduction techniques developed for standard PET can also be applied.An alternative to the proposed delayed window method would be to estimate randoms using singles and coincidence rates, following the theoretical analysis given in Appendix A. However, this approach requires the knowledge of not only the singles rate of 511-keV photons, but also the singles rate of prompt gammas and the coincidence rate of 511-keV photon pairs.
In this paper, we did not perform scatter correction in the SPLIT reconstruction because the objects in the simulation and real experiment are fairly small in size.Both 511-keV photons and prompt gammas can be scattered in PLI.For scattered 511-keV photons, the estimation and correction methods can follow the approaches used in the standard PET image reconstruction for each lifetime-encoded image z T m c by including the estimated scatters in the additive background term r m i in (6).In general, scattered prompt gammas do not affect the reconstrucction of lifetime-encoded images, but they introduce error to the travel time correction.Such effect is already included in the estimated Gaussian blurring function g(τ ) in the pre-reconstruction lifetime fitting.Another problem arises when a prompt gamma is down-scattered into the 511-keV energy window and mis-identified as a 511-keV photon.In this case, it can result in a random tri-coincidence with another unrelated prompt gamma (assuming prompt gamma energy window is sufficiently above 511 keV).The effect of such random events requires further examination.

VII. CONCLUSION
We have developed an efficient positronium lifetime reconstruction method for 3D lifetime image reconstruction using existing TOF PET scanners via time thresholding.The proposed SPLIT method generates a series of lifetime-encoded activity images and estimates the o-Ps lifetime by fitting the lifetime-encoded activity curve for each voxel.We also analyzed the constituents of random events and developed a hybrid strategy to correct for randoms.Both simulation and experimental studies showed that the proposed SPLIT method can reconstruct high-resolution lifetime images with low-variance and outperform the direct TOF-BP method.

APPENDIX A: RANDOMS-TO-TRUES RATIO
Here we perform a theoretical analysis of randoms-totrues ratio in PLI.Let us denote the total activity as A, and the width of the time windows for pairing with another 511-keV photon and a prompt gamma as T 511 and T P G , respectively.For simplicity, a prompt gamma is not considered to be detected in the 511-keV energy window by any chance, which is, however, possible in real cases where a high energy gamma could deposit part of its energy and be identified as a 511-keV event.The discussion is also limited to a relatively small source centered in a cylindrical PET system so that the detection efficiencies can be viewed as unchanged for photons originating from different locations.The detection efficiencies of a single 511-keV photon, a 511-keV photon pair, and a prompt gamma are denoted by P s , P c , and P P G , respectively.In this way, the trues rate is P c P P G A, corresponding to the successful detection of all three photons from a single decay.The randoms rate will be discussed individually for different types of random events (Table I).For type I randoms, the rate of detecting 511-keV pairs is P c A, and the expected number of random prompt-gamma events within the prompt-gamma time window is P P G AT P G .Thus, the rate of type I randoms is P c P P G A 2 T P G .For type I randoms, the rate of detecting one 511-keV event and the prompt gamma associated with it is 2P s P P G A, and the expected number of random 511-keV events within the time window is 2P s AT 511 , where the factor 2 is due to two 511-keV photons per decay.Since the type II.a and II.b randoms have the same rate, the rate of type II randoms is 8P 2 s P P G A 2 T 511 .For type III randoms, the rate of detecting a 511-keV photon is 2P s A, the expected number of random 511-keV events within the time window is 2P s AT 511 , and the expected number of prompt-gamma events within the time window is P P G AT P G .The rate of type I randoms is The above analysis of randoms does not consider the scenario where a prompt gamma is scattered and detected as one (or more) 511-keV photon.To investigate the impact of this situation, we performed two separate simulations of a back-to-back 511-keV source and a prompt-gamma source respectively.The energy spectra are shown in Fig 7 .The ratio of the number of prompt gammas detected in the 511-keV energy window [430,650] keV over the number of true 511-keV photons within the same energy window is 0.087, which may cause error in the randoms-to-trues ratio analysis above.

APPENDIX B: SCANNER SHRINKAGE
The purpose of the scanner shrinkage is to virtually bring the four point sources closer while preserving the spatial and TOF resolutions of the system.Specifically, the detection positions (x, y, z) and TOF difference between two 511 keV photons were first divided by a shrinkage factor η = 30.This shrinkage operation reduced the point sources separation from 16.2 cm to 5.4 mm.However, it also improved the system spatial resolution to 0.21 mm (transaxial) and 1.02 mm (axial), as well as the TOF resolution to 15 ps, much better than those of the original scanner.To restore the spatial and TOF resolution to their original values, the scaled detection positions and TOF difference were then blurred as follows.A Gaussian random variable with a FWHM of 460 ps was added to the scaled TOF value to bring the TOF resolution back to 460 ps.The detection positions were blurred in the axial and tangential directions separately according to the spatial resolution of the original system.
For each threshold T m c , events with lifetime measurements in the window [T 1 , T m c ] are reconstructed by the OSEM algorithm, generating an estimate ẑ(T m c ) of the lifetime-encoded activity concentration z j (T m c ) = x j P T m c ; λ j , A j , m = 1, . . ., M.

Fig. 1 .
Fig. 1.The prompt and delayed time windows.The prompt window is comprised of one 511-keV time window and one prompt-gamma time window.All the boundary values of the time windows displayed are with respect to the time point of the reference 511-keV event.The second 511-keV photon is searched in a range no further than T 511 into the future; the prompt gamma is searched in a range across the reference time point spanned by T a and T b , representing a past and a future time point respectively.T a and T b should satisfy −T a > T M c and −T b < T 1 , respectively.The prompt time window yields true events as well as all the three types of random events.The delayed windows consist of at least one shifted 511-keV or prompt-gamma time window.Each of the delayed time window yields one or two types of random events.

Fig. 2 .
Fig. 2. General workflow of the SPLIT method with random events correction.

Fig. 3 .
Fig. 3.The ground truth lifetime image and the mean and s.d.images of the lifetime images reconstructed by the list-mode ML and the SPLIT methods.

Fig. 4 .
Fig. 4. The lifetime-measurement histograms of the trues, prompts, estimated randoms, and estimated trues (prompts -randoms) in the simulation study.The true events were obtained by examining the eventID output by GATE.

Fig. 5 .
Fig.5.Left: the simulated activity distribution and the activity image of the PLI events reconstructed by standard TOF system matrix.The contrast is not well recovered in the activity image because 44 Sc has a long positron range.Right: the lifetime images reconstructed by the direct TOF-BP and SPLIT methods.The SPLIT method was tested with the (e) true and (f) pre-estimated short lifetime values.The direct TOF-BP method used the true short lifetimes.

Fig. 6 .
Fig. 6.(a) The spatial histogram of the events localized by the direct TOF-BP (dTOF-BP, top row) and activity image reconstructed by OSEM (bottom row) at the original and the virtually shrunk (scale of 1/30) source separations.(b) Reconstructed lifetime images of the four sources by the direct TOF-BP (top row) and SPLIT (bottom row) at different source separations.

2 s
2  P P G A 3 T P G T 511 .Therefore, the overall randoms-to-trues ratio would be randoms rate trues rate= P c P P G A 2 T P G + 8P 2 s P P G A 2 T 511 + 4P 2 s P P G A 3 T P G T 511 P c P P G A = AT P G + P 2 s P c 8AT 511 + 4 A 2 T P G T 511 (A1)which shows that the type I and II randoms-to-trues ratios are proportional to the activity level and the type III randomsto-trues ratio is proportional to the square of activity level.The factor P P c indicates a weak dependence on the scanner sensitivity.

Fig. 7 .
Fig. 7. Energy spectra of a source emitting 511-keV photon pairs and a source emitting 1157-keV (energy of the prompt gamma of Sc-44) photons.

TABLE II COMPARISON
OF THE MEAN ± S.D. OF THE AVERAGE LIFETIME IN THREE ROIS (LEFT AND RIGHT INSERTS, AND THE BACKGROUND) ACROSS 50 REALIZATIONS BETWEEN THE LIST-MODE ML AND SPLIT METHODS

TABLE III RECONSTRUCTED
LIFETIME VALUES IN DIFFERENT ROIS OF THE IMAGES RECONSTRUCTED BY THE DIRECT TOF-BP AND SPLIT METHODS.FOR THE SPLIT METHOD, THE SHORT LIFETIME COMPONENTS WERE FIXED TO EITHER THEIR TRUE VALUES OR THE FITTED VALUES.THE DIRECT TOF-BP METHOD USED THE TRUE SHORT LIFETIMES

TABLE IV MEAN
AND STANDARD DEVIATION OF THE LIFETIME (NS) IN THE FOUR ROIS.THE REFERENCE LIFETIMES WERE OBTAINED BY FITTING THE LIFETIME-MEASUREMENT HISTOGRAM OF THE EVENTS ORIGINATED IN THE NEIGHBORHOOD OF EACH SOURCE