A Conjugate Bayesian Framework for Fast 3D Positronium Lifetime Estimation with a Partial System Matrix

Imagine you are trying to figure out how long a specific type of firefly lives inside a giant, dark, 3D maze. You can't see the fireflies directly, but you have thousands of sensors on the walls of the maze. Every time a firefly dies, it sends out a flash of light (a photon) that hits the sensors.

Positronium Lifetime Imaging (PLI) is like trying to map exactly how long these "fireflies" (which are actually subatomic particles called positronium) live in every single tiny cube of space inside a patient's body. This is super useful because the lifespan of these particles changes depending on the "environment" they are in (like whether they are in healthy tissue or a tumor).

However, doing this in 3D is a massive headache for computers. Here is why, and how this paper solves it.

The Problem: The "Too Many Possibilities" Trap

Imagine the maze has 1,000 sensors. If you try to calculate every possible path a firefly could take from every point in the maze to every possible sensor, the number of calculations is astronomical. It's like trying to write down every possible sentence in the English language just to find the one sentence you actually heard.

In the real world, we only catch a tiny fraction of these "firefly deaths." Most of the possible paths never happen. But traditional computer methods try to calculate the math for all the impossible paths anyway, which crashes the computer or takes days to finish.

The Solution: The "Smart Filter" Framework

The authors of this paper built a new, super-fast way to do this math. They used three clever tricks, which we can explain with analogies:

1. The "Only Look at What Happened" Rule (Partial System Matrix)

Instead of trying to calculate the math for every possible sensor combination in the universe, their new method says: "Let's only look at the specific paths that actually happened."

Analogy: Imagine you are a detective trying to solve a crime. Instead of interviewing every person in the city (most of whom were nowhere near the crime), you only interview the people who were actually seen at the scene.
Result: This cuts the amount of data the computer has to process by a huge amount, making it possible to run the calculation in seconds instead of days.

2. The "Fractional Credit" System (Event-to-Voxel Weighting)

When a sensor sees a flash, it doesn't know exactly which tiny cube (voxel) in the 3D maze the firefly came from. It's like hearing a noise in a crowded room and not knowing exactly which person made it.

Analogy: Instead of guessing "It was definitely Person A" or "It was definitely Person B," the computer says, "There is a 70% chance it was Person A and a 30% chance it was Person B." It gives "fractional credit" to everyone who might have made the noise, based on where they were standing and how loud the noise was.
Result: This allows the computer to use every single piece of data it has, even when it's not 100% sure where it came from, without getting confused.

3. The "Instant Calculator" (Conjugate Bayesian Update)

Usually, to figure out the average lifespan of the fireflies, computers have to guess, check, guess again, and check again (iterative loops) until they get the answer right. This is slow.

Analogy: Imagine trying to find the average height of a group of people.
- Old Way: You guess a height, measure everyone, realize you were wrong, guess again, measure again... repeat 10 times.
- New Way: The authors found a special mathematical "shortcut" (a conjugate update) that lets them calculate the answer in one single step using a formula, just like using a calculator instead of counting on your fingers.
Result: This turns a process that took over a minute into one that takes a fraction of a second.

What Did They Prove?

They tested this new method in two ways:

Computer Simulation: They created a fake 3D world with known firefly lifespans. Their new method found the correct answer in 2.76 seconds, while the old, standard method took 74 seconds on the same computer.
Real Hardware: They used a real prototype scanner (J-PET) and a test phantom (a fake body made of plastic spheres). They successfully mapped the lifespans of over 230,000 tiny cubes in just 3 seconds.

The Bottom Line

This paper is like upgrading from a manual typewriter to a word processor with "Auto-Correct" and "Search."

They didn't invent a new type of firefly or a new sensor. Instead, they invented a smarter way to organize the data. By ignoring the impossible paths, sharing the credit for uncertain events, and using a mathematical shortcut to solve the final equation, they made it possible to create detailed 3D maps of how particles behave inside the body in real-time. This could eventually help doctors spot diseases like cancer much faster and more accurately.

1. Problem Statement

Positronium Lifetime Imaging (PLI) extends conventional Positron Emission Tomography (PET) by utilizing the time interval between positron emission and annihilation as a contrast mechanism. This provides information about the tissue microenvironment. However, performing voxel-wise lifetime estimation in fully 3D settings presents significant computational challenges:

Data Sparsity vs. System Complexity: The number of feasible detector-time channels (combinations of detector pairs and time-of-flight bins) grows rapidly ( $O(N_{sensor}^2 \times N_{TOF})$ ), creating a massive system matrix.
Observation Limitation: In practice, only a tiny fraction of these theoretical channels are actually observed in the data (sparse triple coincidences).
Computational Cost: Traditional iterative optimization methods (like Maximum Likelihood Estimation) applied to the full system matrix or even large sparse matrices are computationally prohibitive for large 3D volumes, often requiring hours to converge.

The paper addresses the need for a scalable, computationally efficient framework to estimate 3D positronium effective lifetimes from sparse triple-coincidence data.

2. Methodology

The authors propose a two-stage statistical estimation framework that replaces the full system matrix with a "partial" one and utilizes a conjugate Bayesian update.

A. Partial System Matrix (Forward Model)

Instead of modeling all possible detector-time channels, the forward model is restricted strictly to the observed channels ( $C^+$ ) present in the acquired data.

Definition: The system matrix $H_{c,j}$ is defined as the conditional probability $P(C=c | J=j, C \in C^+)$ , where $c$ is an observed channel and $j$ is a voxel.
Normalization: The matrix is normalized over the observed channels for each voxel ( $\sum_{c \in C^+} H_{c,j} = 1$ ).
Benefit: This drastically reduces memory usage and computational cost while preserving the Poisson data model for the retained triple coincidences. The model estimates the distribution of retained detected events, not the total physical decays.

B. Two-Stage Estimation Strategy

Stage 1: Activity Reconstruction
- Estimate the expected number of retained triple coincidences per voxel ( $\tilde{f}_j$ ).
- Method: Expectation-Maximization (EM) algorithm maximizing the Poisson likelihood.
- Advantage: Because the system matrix is normalized over observed channels, the sensitivity term in the EM update is unity, simplifying the backprojection step.
Stage 2: Lifetime Estimation
- Conditional on the activity estimate $\tilde{f}$ , estimate the voxel-wise effective annihilation rate ( $\lambda_j$ ) using the list-mode lifetime data ( $\tau_k$ ).
- Event-to-Voxel Weighting: Since the origin voxel of an event is unknown, the method calculates posterior weights $\pi_{k,j}$ representing the probability that event $k$ originated from voxel $j$ , given the observed channel and the reconstructed activity.
- Two Estimators Compared:
  1. Maximum Likelihood (L-BFGS-B): Uses the weighted events to maximize the log-likelihood iteratively.
  2. Conjugate Bayesian Estimator (Proposed):
    - Assumes a Gamma prior for $\lambda_j$ .
    - Replaces unknown latent indicators with the calculated posterior weights $\pi_{k,j}$ .
    - Computes sufficient statistics (weighted counts and weighted time sums) in a single pass.
    - Derives a closed-form posterior mean for the rate:
      $\hat{\lambda}_j = \frac{\alpha_0 + \sum_k \pi_{k,j}}{\beta_0 + \sum_k \pi_{k,j}\tau_k}$
    - This avoids iterative optimization entirely.

3. Key Contributions

TOF-Aware Partial System Matrix: A formulation that restricts the forward model to observed detector-time channels, making 3D reconstruction computationally tractable for sparse data while maintaining statistical rigor.
Principled Event Attribution: A weighting scheme that treats the source voxel as a latent variable, allowing for fractional attribution of events across voxels based on detector geometry and reconstructed activity.
Conjugate Bayesian Surrogate: A closed-form analytic estimator for voxel-wise effective rates. It provides a fast, stable alternative to iterative optimization, offering uncertainty quantification (via the posterior distribution) without the computational burden of L-BFGS-B.

4. Results

The framework was evaluated using both simulated data and experimental data from a J-PET prototype scanner with a NEMA phantom.

Simulation Study (4,056 voxels):
- Accuracy: The Bayesian estimator accurately recovered the effective-rate map, matching the spatial patterns and dynamic range of the ground truth.
- Speed: The analytic Bayesian estimator required 2.76 seconds, compared to 74.46 seconds for 10 iterations of L-BFGS-B on the same CPU (approx. 27x speedup).
- Error: Both methods showed median relative errors near zero. The Bayesian method exhibited tighter error distributions (lower variance).
Experimental Study (J-PET, 234,375 voxels):
- Scale: Successfully estimated a lifetime map for a large volume (234k voxels) from ~364,000 retained events.
- Performance: The full estimation (activity + lifetime) completed in approximately 3 seconds on a single CPU.
- Qualitative Results: The reconstructed activity map correctly identified the $^{44}$ Sc-filled spheres. The lifetime map provided a stable single-component surrogate summary of the timing signal, consistent with expected effective lifetimes (~1.3 ns), though interpreted as a surrogate due to the simplified single-component model.

5. Significance and Conclusion

This work establishes a practical computational foundation for large-scale 3D Positronium Lifetime Imaging.

Scalability: By restricting the system matrix to observed channels, the method overcomes the memory and speed bottlenecks that previously hindered 3D PLI.
Efficiency: The conjugate Bayesian update eliminates the need for slow iterative optimization, enabling near real-time 3D reconstruction.
Robustness: The framework provides stable estimates even in sparse data regimes where individual voxels may have few detected events, thanks to the probabilistic event weighting.
Future Outlook: While the current study uses a simplified single-component exponential model, the framework is designed to be a scalable surrogate. It paves the way for future integration of more complex multi-component decay models and detector response corrections in clinical and research settings.

In summary, the paper demonstrates that combining a partial system matrix with a conjugate Bayesian update makes fully 3D positronium lifetime imaging computationally feasible, accurate, and fast enough for practical application.