A dynamic mPET scan records the spatiotemporal distribution of a mixed uptake of N tracers within a living organism. The multi-tracer activity concentrations (i.e., multi-tracer TAC¹) of an image voxel (or in a ROI) at time t can be modeled as the linear superposition of the pharmacokinetic model of each tracer (16, 17, 19)(1)CMulti(t;κ)=VBS(t)+(1−VB)∑n=1N[A(n)(t)⊗R(n)(t;κ(n))]e−λ(n)t,where κ is a vector that contains all tracer kinetic parameters, R(n)(t) is the impulse response function of the nth tracer, λ(n) denotes the rate of radioactive decay, A(n)(t) is the tracer concentration in plasma (i.e., AIF), S(t) is the total activity concentration in whole blood, VB∈(0,1] is the fractional volume of blood in the tissue, and ⊗ denotes the convolution operator. The mPET image intensity at voxel j in time frame k is then given by Equation 2:(2)xkMulti(κj)=∫tk,stk,eCMulti(τ;κj)dτ,where tk,s and tk,e represent the start and end points of frame k. The expectation of the projection mPET data y¯k(κ) in time frame k with respect to the dynamic mPET image xk(κ) can be expressed by Equation 3:(3)y¯kMulti(κ)=PxkMulti(κ)+rk,where the (i,j)th element of the system matrix P∈RI×J is the probability of detecting an event originating in voxel j by detector pair i, I and J are the total number of detector pairs and image voxels, respectively, and rk is the expectation of scattered and random events in the kth frame.

The overview of the proposed deep learned post-separation framework for triple-tracer PET image separation is illustrated in Figure 1. In our proposed DL-based method, the triple-tracer activity images (voxel-wise TACs) were used as the network input and the single-tracer activity images were used as the training labels.

We employed an architecture based on the CED, which has been widely used in the DL-based mPET image separation (31, 33, 38–40). In the proposed network, the encoder branches consist of the repeated application of two 3×3 2D convolution layers, each followed by a batch normalisation (BN) and a parametric rectified linear unit (PReLU), in addition to a max-pooling layer for downsampling, followed by the BN and the PReLU. Each of the decoder branches consists of a 3×3 2D transposed convolution layer for upsampling and two 3×3 2D convolution layers, each followed by the BN and the PReLU, and a 1×1 2D convolution layer at the end of each decoder branch. In addition, we activated the output layer using a ReLU to enforce the non-negativity constraint on the separated single-tracer activity images. The number of trainable parameters for the proposed CED 2D is approximately 2.1×106. The mean squared error (MSE) loss is applied to the activity images of each tracer, and their sum Ltotal is used as the loss function for network training, which is given by Equation 4:(4)Ltotal=∑n=13∑s=1S1S‖C∘s(n)−C^s(n)‖22,where C∘(n) denotes the activity images of the nth tracer in the network output, C^(n) is the activity images of the label single-tracer, and S is the total number of training pairs.

We have focused on the simulation study of the triple-tracer myocardium PET image separation based on the combination of [F18]FDG, Rb82, and [Tc99m]sestamibi with decay constants λ[F18]=log⁡(2)/109.7min−1, λ[Rb82]=log⁡(2)/1.26min−1 and λ[Tc99m]=log⁡(2)/52min−1. A simulation was performed to assess the performance of the proposed DL-based triple-tracer separation method. A myocardium dataset was obtained from the 3D MRI cardiac scans provided in the M&Ms challenge² (43). Myocardium segments were extracted from two non-continuous slices (short-axis view), which were selected from one of the frames of each 3D dynamic MRI image. A total of 140 myocardium segment images were obtained. Each myocardium segment image was resized to 128×128 with a voxel size of 2.602×2.602 mm2, and further randomly divided into 4 to 14 sub-regions with well-defined boundaries. The ground-truth kinetic parameter for a given sub-region was sampled from a Gaussian distribution with mean values derived from the literature (44–46) (see Table 1)³ and coefficient of variation equal to 0.1 to simulate heterogeneous variation within the whole myocardium region (absolute values were taken after sampling). A simulated ground-truth K1 parametric map of FDG is shown in Figure 2A (for plotting only, we restricted the voxel values to [0.4, 0.7] to demonstrate the sub-region segments), along with the pre-defined myocardium ROI. The AIFs of each tracer were generated along the shape of the AIFs from the literature (44–46) using Feng’s input function model (47). To further simulate the population variation in the myocardium dataset, the parameters of Feng's input function model were also modelled as a Gaussian variable with a coefficient of variation equal to 0.1 (39), and the absolute values were taken after sampling.

The ground-truth single-tracer voxel-wise TACs were generated from the simulated parametric maps using the irreversible two-tissue compartment model for FDG and Rb, and the reversible two-tissue compartment model for sestamibi. The single-tracer voxel-wise TACs were then summed up together to form the ground-truth triple-tracer voxel-wise TACs. In this study, we propose a realistic protocol for the dynamic triple-tracer PET scan, which was conducted for 60 min after the [F18]FDG injection (with an initial 60 min off-bed time). Rb82 was injected 5 min after the start of the dynamic PET scan, followed by the injection of [Tc99m]sestamibi with a 10 min delay. An example of the simulated AIFs and the ROI TACs without decay correction based on the proposed scan protocol is shown in Figure 3.

For the reconstruction of the simulated data, we modelled a GE Discovery ST PET-CT scanner with a system sensitivity of ∼2 cps/kBq in 2D mode (48). Noise-free sinogram data were generated by forward-projecting the dynamic ground-truth images (for both triple-tracer and single-tracer) using a pre-calculated system matrix. A 20% uniform background was included as a simple model for the mean of the random and scatter background events. A scaling factor was applied to this projected data in order to generate the mean count levels in the sinogram, prior to the introduction of Poisson noise into each sinogram bin. This scaling factor was chosen to obtain datasets containing mean total counts in each 2D sinogram, where the mean total counts were modelled based on the system sensitivity of the scanner (Figure 2B).

The dynamic PET scan (60 min for triple-tracer and single-tracer) was divided into 28 time frames: 1×5 min, 4×0.25 min, 2×0.5 min, 3×1 min, 1×2 min, 1×3 min, 4×0.25 min, 2×0.5 min, 3×1 min, 1×2 min, 1×3 min, 3×5 min, 2×10 min. Dynamic images were reconstructed using the MLEM algorithm (initialised by uniform images) with 128 iterations without post-smoothing. The reconstructed images were frame-length corrected and thus equivalent to reconstructed tracer activity images (voxel-wise TACs).

For the proposed DL-based method, the dynamic triple-tracer noisy MLEM reconstruction was used as the network input, and the dynamic single-tracer noisy MLEM reconstructions were used as the training labels. We compared the separation performance of the triple-tracer (TT) with the dual-tracer (DT) and the single-tracer (ST) using the proposed DL-based method. For the dual-tracer separation, dynamic dual-tracer ([F18]FDG and Rb82) noisy MLEM reconstructions were used as the network input and the loss function Ltotal (with n=2) was applied to [F18]FDG and Rb82 only, and dynamic noisy single-tracer MLEM reconstructions were used as training labels. For the single-tracer case, the dynamic noisy single-tracer MLEM reconstruction was used as the network input. To have a fair comparison, the dynamic noisy single-tracer MLEM reconstructions were also used as training labels (with n=1 in Ltotal). All network training was performed in the same manner. In total, 120 simulated data examples were used for training, 10 for validation and 10 for testing. The network parameters were initialised using the Xavier initialisation. The Adam algorithm (49) was used with a learning rate of 5×10−4 and a batch size equal to 8 for network training. To prevent overfitting, all networks were trained for 1,500 epochs with early stropping if there was no improvement in the validation metrics. The network training and evaluation steps were implemented in PyTorch, on a PC with an NVIDIA GeForce RTX 3,090 GPU.

We also compared the DL-based method with the voxel-wise MTCM-based method (v-MTCM) (14, 16). The v-MTCM method estimates the single-tracer TACs by fitting the triple-tracer kinetic model as shown in Equation 1 (with known AIFs of each tracer⁴) to the measured triple-tracer TACs, i.e., the dynamic triple-tracer noisy MLEM, using voxel-wise weighted least squares (VWLS). Note that the fitting was only performed on the frames starting at 60 min after FDG injection, i.e., on the dynamic triple-tracer images obtained within the 60 min triple-tracer PET scan. The time frame durations were used as weighting factors to compensate for non-uniform temporal sampling (16). The trust-region-reflective algorithm was used to perform the VWLS fitting. The stopping criteria were set such that the optimisation procedure was stopped when the relative termination tolerance of the objective function was less than 1×10−8 or the maximum number of iterations (1,600 iterates) was achieved (21). The initial values of the kinetic parameters κ were set to be 0.01 for all voxels and the values of the lower bounds for each parameter were set to be 1×10−5, while the values of the upper bounds for VB, K1 and k2−k4 were set to be [1,5,2,1,1,1], respectively.

Both the DL-based and the MTCM-based methods can be implemented at the ROI level for triple-tracer ROI-TAC separation. The triple-tracer ROI TACs were extracted from the dynamic triple-tracer noisy MLEM followed by a 1D TAC separation. For the DL-based method, we simply replaced the 2D modules in the proposed network shown in Figure 1 with their 1D versions (CED 1D). The number of trainable parameters for the CED 1D is approximately 7.7×105. The triple-tracer ROI TACs were used as the network input and the single-tracer ROI TACs extracted from the dynamic single-tracer noisy MLEM were used as the training labels. To perform the MTCM-based method at the ROI level (ROI-MTCM), the triple-tracer compartment model was fitted to the extracted triple-tracer ROI TACs to recover the ROI TACs of each tracer using the single-voxel WLS. Both the CED 1D and the ROI-MTCM methods were implemented in the same manner as described for the voxel-level separation.

Figure 4 shows the single-tracer noise-free MLEM and the separated single-tracer activity images by different separation methods for frame 16 (a 15-s frame, at 1 min after the injection of [Tc99m]sestamibi). It can be observed that the triple-tracer separation using the CED 2D (with 120 training examples) can substantially reduce the image noise compared to the triple-tracer v-MTCM separation. Figure 5 shows the quantification results of the separated single-tracer activity images (all time frames were considered). The performance of the triple-tracer separations by using different methods is shown in the top row. For the single-tracer activity images, lower bias and SD were found by fitting the compartment model to the single-tracer noisy MLEM voxel by voxel (v-STCM) for FDG, Rb and sestamibi (moving from the yellow cluster to the black cluster), which is consistent with the visual impression of the activity images shown in Figure 4 (by comparing the second and third columns). Although the single-tracer CED 2D further reduces the SD compared to the v-STCM, a higher bias was still obtained due to the small number of training pairs and the prior of the proposed network.

For the considered triple-tracer separation, the situation worsens. High bias was obtained for each tracer using the v-MTCM method for the noisy triple-tracer separation compared to the single-tracer cases. The main reason is that the MTCM-based method is sensitive to noise and may fall into local minima even when the noise level is low (21). This can be verified by comparing the performance of the v-MTCM method for noise-free triple-tracer separation,⁵ where the SD is much lower, but the bias level is still higher than the single-tracer noisy v-MTCM. The DL-based method, CED 2D, sufficiently reduces the SD and achieved lower bias compared to the v-MTCM method for noisy triple-tracer separation even though the network was trained using noisy labels. The CED 2D uses both spatial and temporal information for triple-tracer separation, while the v-MTCM method only uses temporal information and is highly dependent on prior information, such as the time delay intervals between tracer injections (33, 42). More importantly, the CED 2D with MSE loss learns to output the mean of all plausible noisy explanations when it is trained using noisy labels, and thus the proposed network implicitly learns to denoise the output images (33). However, the CED 2D for triple-tracer separation fails to reach the same level of bias and SD as the CED 2D for single-tracer using the same number of training examples, indicating that the triple-tracer separation task is still challenging for the proposed DL-based method.

Dual-tracer ([F18]FDG and Rb82) separation was also investigated. The bias and SD trade-off is shown in the bottom row of Figure 5. The CED 2D results in lower SD but higher bias compared to the v-MTCM for noisy dual-tracer separation for both FDG and Rb. The lower bias level in the v-MTCM separation is due to several factors: (i) the very short half-life of Rb82 (76-s), (ii) the injection of Rb82 was delayed by 5 min after the start of the dynamic PET scan, i.e., only the [F18]FDG signal was measured in the first 5 min, and (iii) the activity concentration of [F18]FDG tends to become stable (see the rather flat red line in Figure 3). These factors also make the dual-tracer signals become much easier to disentangle compared to the triple-tracer separation task, leading to lower bias and SD in the dual-tracer separation compared to the triple-tracer case for both the v-MTCM and the CED 2D. We also noted that the bias and SD of Rb are higher than those of FDG and sestamibi for the single-tracer, the dual-tracer, and the triple-tracer cases. This is because the MLEM reconstructions of Rb were extremely noisy due to the low-count level caused by the short half-life of Rb82 and the short time frame duration in the early time frames.

We retrained each CED 2D using different sample sizes to assess the impact of the number of training examples on triple-tracer separation. Figure 6 shows that the NRMSE values of the separated single-tracer activity images decrease as the increasing number of training examples (from 8 to 120). With 120 training examples, the CED 2D triple-tracer separation achieved ∼12%, ∼15% and ∼12% NRMSE values for FDG, Rb and sestamibi, respectively. However, these errors are still higher than those of the CED 2D for single-tracer using the same number of training examples (120 pairs), which again shows the bottleneck of the proposed DL-based separation method.

Parametric imaging was also performed on the separated activity images. The v-STCM method was used to estimate the parametric maps from the separated single-tracer images obtained from the CED 2D. Note that in the v-MTCM method, the parametric images of each tracer were separated before recovering the single-tracer activity images, eliminating the need for post-estimation. The parametric maps recovered from the single-tracer noise-free MLEM⁶ were used as reference images. Figure 7 shows the parametric images of [F18]FDG for a static image corresponding to the last 20-min interval obtained by frame integration, Rb82 delivery rate k3, and [Tc99m]sestamibi delivery rate K1. Compared with the parametric images estimated from the triple-tracer v-MTCM separated activity images, the CED 2D results are more similar to the reference images.

Figure 8 shows the corresponding bias and SD trade-off of the separated parametric images. Without the impact of noise, the parametric maps obtained from the noise-free v-MTCM separation achieved much lower bias and SD compared to the noisy v-MTCM separation for both the triple and dual-tracer separations. The parametric maps obtained from the CED 2D separation exhibit a dramatic reduction in SD compared to the noisy v-MTCM separation. In addition, lower bias and SD of the parametric maps were achieved in the dual-tracer separations compared to the triple-tracer separations for both the v-MTCM and the CED 2D (see Subsection 4.1 for detailed discussion). However, the bias of the static FDG and sestamibi K1 estimated from the triple-tracer and dual-tracer separations using the v-MTCM method is lower than that of the CED 2D, indicating the lack of training examples and the weak inductive prior of the CED 2D. The Rb k3 images generally have higher bias and SD compared to the FDG static and sestamibi K1 images using the v-MTCM, demonstrating that the separation and estimation of Rb k3 is still challenging. However, a clear reduction in bias and SD for k3 was achieved by using the CED 2D, which is consistent with the observations in Figure 7.

The separated ROI TACs using the MTCM-based and DL-based methods were also assessed. Figure 9 shows the mean separation results of a single test example over 20 different noise realisations. The reference single-tracer ROI TACs extracted from the single-tracer noise-free MLEM (dashed lines), the separated ROI TACs using the voxel-level methods (v-MTCM and CED 2D) and the ROI-level methods (ROI-MTCM and CED 1D) are presented. Both the CED 2D and the CED 1D generally resemble the reference TACs of each tracer compared to the v-MTCM and the ROI-MTCM.

The NRMSE of the separated TACs are shown in Figure 10. For the MTCM-based separation methods, the v-MTCM results in much lower NRMSE for all tracers compared to the ROI-MTCM, indicating that the separation performance was enhanced for the MTCM-based method by considering the voxel-wise separation in the myocardium ROI. On average, a significant reduction in NRMSE was achieved by the CED 2D and the CED 1D compared to the v-MTCM and the ROI-MTCM for all tracers. The NRMSE of the CED 1D is higher than that of the CED 2D for FDG, while it is lower for Rb and sestamibi. We observed that the ROI TAC of FDG became stable (see the rather flat red dashed lines in Figure 9). The CED 2D implicitly learned to denoise the output images, which is beneficial for the separation of the rather flat FDG ROI TAC. However, the Rb and sestamibi ROI TACs each contain a peak at the very beginning, which is more difficult to recover in the separation compared to the FDG ROI TAC. In this case, the CED 1D results in a better separation performance. The CED 1D learns a direct mapping to separate single-tracer TACs from triple-tracer TACs in the time domain, which is more effective, whereas the CED 2D is an indirect method where single-tracer TACs are extracted from single-tracer images after the separation in the image domain. However, the NRMSE of the separated TACs using either CED 1D or CED 2D in general fails to reach the same level as the single-tracer noisy MLEM as highlighted by the dashed yellow lines.