Design of a generic method for single dual-tracer PET imaging acquisition in clinical routine

Using different tracers in positron emission tomography (PET) imaging can bring complementary information on tumor heterogeneities. Ideally, PET images of different tracers should be acquired simultaneously to avoid the bias induced by movement and physiological changes between sequential acquisitions. Previous studies have demonstrated the feasibility of recovering separated PET signals or parameters of two or more tracers injected (quasi-)simultaneously in a single acquisition. In this study, a generic framework in the context of dual-tracer PET acquisition is proposed where no strong kinetic assumptions nor specific tuning of parameters are required. The performances of the framework were assessed through simulations involving the combination of [18F]FCH and [18F]FDG injections, two protocols (90 and 60 min acquisition durations) and various activity ratios between the two injections. Preclinical experiments with the same radiotracers were also conducted. Results demonstrate the ability of the method both to extract separated arterial input functions (AIF) from noisy image-derived input function and to separate the dynamic signals and further estimate kinetic parameters. The compromise between bias and variance associated with the estimation of net influx rate K i shows that it is preferable to use the second injected radiotracer with twice the activity of the first for both 90 min [18F]FCH+[18F]FDG and 60 min [18F]FDG+[18F]FCH protocols. In these optimal settings, the weighted mean-squared-error of the estimated AIF was always less than 7%. The K i bias was similar to the one of single-tracer acquisitions; below 5%. Compared to single-tracer results, the variance of K i was twice more for 90 min dual-tracer scenario and four times more for the 60 min scenario. The generic design of the method makes it easy to use for other pairs of radiotracers and even for more than two tracers. The absence of strong kinetic assumptions and tuning parameters makes it suitable for a possible use in clinical routine.


Introduction
Personalized medicine is one of the significant aims of oncology. Thanks to recent technological breakthroughs, the fine characterization of tumors allows for the identification of diagnostic biomarkers which can be useful in the context of prognosis or therapeutic response. Positron emission tomography (PET) is a powerful tool that allows characterizing the tumor and its microenvironment in an in vivo and non-invasive way. The use of various PET radiotracers can provide complementary information at the voxel level about intra-tumoral heterogeneities due to the different biological processes of radiotracers (Bailly et al 2019).
For instance, [ 18 F]FDG is a glucose analog and allows to examine glucose metabolism while [ 18 F]FCH is related to the synthesis of plasma membranes and then linked to cells proliferation. It was demonstrated that information provided by [ 18 F]FCH and [ 18 F]FDG complement each other in different clinical indications such as multiple myeloma (Mesguich et al 2020). In another study, it was shown that PET imaging using [ 18 F]FCH and [ 18 F]FDG provided more accurate detection of low-grade glioma than using only one of these radiotracers (Fu et al 2016). Therefore, lesions may be better characterized through PET imaging using two or more radiotracers for the same subject.
Usually, imaging the distribution of two different radiotracers requires performing two sequential singletracer PET acquisitions with a time delay between injections allowing for the complete decay of the first tracer. However, internal and external movements impede the accurate registration of the two PET images. In addition, physiological changes may occur between the two acquisitions.
To address these issues, a single PET acquisition with (almost) simultaneous injections of two radiotracers has been proposed (Huang et al 1982, Koeppe et al 2001, Black et al 2009, Andreyev and Celler 2011, Kadrmas et al 2013. The subsequent challenge of this dual-tracer PET technique is to extract the signal corresponding to each of the injections from the PET kinetics resulting from the co-injection, in order to derive the biological parameters specific to the radiotracer. In the case of radiotracers with different isotopes (e.g. [ 18 F] and [ 11 C]), Huang et al (1982) took advantage of the difference in decay rates. Isotopes emitting additional high-energy simultaneously with a positron can also be used to discriminate each tracer (Andreyev and Celler 2011). The signal separation can be achieved by fitting the tracers pharmacokinetics using compartmental models. In that respect, Koeppe et al (2001) (Kadrmas et al 2013). Bell et al (2017) proposed a generic Spectral Analysis (SA) dynamic model (Cunningham and Jones 1993) with an arbitrarily high number of spectral functions and an ℓ 1 regularization term to promote sparsity in the selection of contributing spectral functions. In this optimization scheme, the strength of the sparsity constraint has to be tuned. This approach requires the measurement of the arterial input function (AIF), which generally relies on arterial blood sampling.
Bell and colleagues considered this latter to be known without an explicit method to derive it. In contrast, Kadrmas et al (2013) proposed a chemical separation of input functions (IFs) from the arterial samples. An alternative to arterial blood sampling, which is invasive, is to derive the input function from the image itself, usually using a region of interest placed on arteries (Hu et al 2020). This image-derived input function (IDIF) is efficient but requires adjustment and modeling of the AIF.
With the goal of performing dual-tracer PET imaging in clinical routine, we developed a generic simple and robust framework aiming at computing pharmacokinetics parameters of two distinct co-injected radiotracers, with minimal assumptions.
In the fitting process, AIFs were extracted and separated from the noisy IDIF signal using a model presented in section 2.1. Time activity curves (TAC)s were then recovered by coupling the SA model together with the nonnegative least-square (NNLS) algorithm as detailed in section 2.2. Simulations of TACs with [ 18 F]FDG and [ 18 F]FCH were used to extensively characterize the performance of the proposed method. Eventually, preclinical experiments were conducted to consolidate the results obtained from the simulations.

Methods
All Python programming scripts implementing the methods presented hereafter are publicly available at https://gitlab.com/ntaheri/dual-tracer-codes.
2.1. Input function 2.1.1. AIF fitting In dual-tracer PET, arterial blood activity is a mixture of each tracer activity. Therefore, since AIF is required for tracer kinetics modeling, it has to be independently determined for each tracer from the mixed signal. For a single injection, the AIF can be modeled by a linear function from the injection time to the peak concentration time, and by a sum of three exponential terms after the peak (Vriens et al 2009). For a dual injection, considering t p 1,2 and τ 1,2 the peak times and the injection times of respectively the first and second tracer, the dual-tracer AIF can be modeled with the following equations: Calculating the exponential parameters A 1 i , A 2 i and also g 1 i , g 2 i , for i = 0, 1, 2, allows to reconstruct the AIF for the first and second tracer; i.e. AIF Tracer1 (t) and AIF Tracer2 (t). In addition, it is assumed that τ 1 = 0.
2.1.2. AIF fitting from IDIF curve As an alternative to invasive AIF, IDIF could be determined by measuring the activity concentration taken in the large arteries or heart on dynamic PET images. Dynamic PET data consists of a series of images being the average of the tracer concentrations over time frames. So modeling the AIF from the IDIF consists in computing a series of integrals over time frames. It can be expressed as the following equation for a specific frame n and a time length l[n]:  with the respective radioisotope half-life T 1 2 x .
2.2. PET kinetics 2.2.1. Modeling the tissue TACs A TAC of a tracer in a tissue can be described by a compartment model. This model is governed by a set of differential equations involving the concentration of the tracer in each compartment and exchange rates between compartments. These equations can be used either to simulate the tissue TACs or to fit the model parameters from measured dynamic PET data. The compartment model needs some a priori assumptions on the tracer kinetics behavior in tissue (e.g. number of compartments, reversibility). In this study, the two tissue compartment model (2TCM) was only used to simulate the noise-free tissue TACs of irreversible kinetics (the equations can be found in appendix A.5). In this work, the PET TACs were determined using SA since it does not require any assumptions on tracer kinetics, contrary to the compartment model.

SA with NNLS algorithm
To build a generic framework without requiring any a priori information about tracer kinetics, we used the SA to model the dynamic PET (Cunningham and Jones 1993). Exploiting the SA and arguments from Slawski and Hein (2013), together with the NNLS algorithm for the fitting process, allowed not only to dispose of any tuning parameters but also to retain the sparsity of the spectral coefficients.
The tissue concentration of a dual-tracer at time t, C(t) = C Tracer1 (t) + C Tracer2 (t − τ 2 ) is modeled as a sum of N possible tissue responses. This model is a linear combination of a convolution of AIF(t), by a set of exponential terms such that: The values of θ are selected to cover an appropriate spectral range for both smooth and transient events of the tissue tracer concentration. For instance, θ = 0 is considered as the lowest frequency component and is related to the slowest kinetics. On the other hand, very large values of θ are considered as high-frequency components closer to the AIF. The convolution operator is denoted by ⊗. The values α i and θ i are assumed to be non-negative due to the constraint of first-order tracer kinetics. The dynamic responses (basis functions) are defined as (Bell et al 2017). After setting N and θ, the values of α are estimated from the AIF and tissue TAC by a weighted NNLS algorithm. The estimated vector α is sparse as a few numbers of the coefficient are non-negative.
A timing resolution of 0.1 s was used to build the input function and compute the basis functions. The basis functions included the input function (ψ N ), the integral of the input function (ψ 0 ) and N = 40 functions ψ using a set of θ i as described in Veronese et al (2010). The θ values are logarithmically spaced within the range of , duration of the first frame). Thanks to the properties of NNLS, higher values of N are not expected to change the results but will suffer from higher computational burden. Consequently, the value of N is selected to give the best results with less computation.

Patlak modeling
In this work, we especially looked at the net influx rate K i . As such, a Patlak modeling (Patlak et al 1983) using the TACs of the tracers estimated from the SA coefficients was also conducted with the hypothesis that it would help to get a more robust estimation of derived macroparameters.

Simulations
The IDIF measurements have been simulated using the parameters extracted from the real single-tracer acquisitions described in section 2.4. The obtained IDIF parameters for [ 18 F]FCH are t p = 3.66s (time of peak), The IDIF data were generated by giving the mentioned parameters to model (2). An example of noisy IDIF, noise-free AIF and estimated AIF signals are displayed in figure 1 (Schaefferkoetter et al 2017).
For the case of [ 18 F]FDG, it is recommended to obtain the kinetic parameters from dynamic acquisition 30-60 min after tracer administration (Hu et al 2020). Therefore, this scenario was considered favorable. To study the quality of a shortened scan, which would be of interest for patient comfort, 60 min acquisitions were also simulated with reversed order of injections ( The timing frames of single-tracer simulations were [2 × 10 s, 2 × 20 s, 4 × 30 s, 3 × 60 s, 2 × 120 s, 5 × 240 s, 1 × 600 s, 1 × 900 s] for the 60 min simulations and [2 × 10 s, 2 × 20 s, 4 × 30 s, 3 × 60 s, 2 × 120 s, 5 × 240 s, 2 × 450 s, 3 × 900 s] for the 90 min simulations. For dual-tracer simulations, this framing was repeated at the injection time of the second tracer.
Gaussian noise was added to the simulated TACs. The weight w of frame n was calculated using the radioactive decay constant λ, the start time of the frame T s and the duration of the frame T d as: The standard deviation of the Gaussian distribution σ of frame n was calculated using the weighted ratio of noise-free signal s[n] multiplied by the value ν (noise intensity) as follows: In order to have realistic measurements, ν = 0.35 was selected. The impact of different noise levels is presented in appendix A.3. To generate noisy TACs, the noise-free TAC simulations were contaminated as follows: TAC noisy = TAC noise−free + σ. One thousand realizations of noisy TACs were generated for each experiment.
To study the impact of the activity of the radiotracers on the estimation of macroparameters, different ratios of activity (initial activity of the first tracer over the initial activity of the second tracer, corrected for decay) were considered. The studied ratios were 0.25, 0.5, 1, 2, 4, considering the clinical constraint that the summation of the activities is equal to two. For single-tracer measurements, one quantity of the radiotracer was applied.
The AIF estimation derived from noisy IDIF was compared to the noise-free AIF using the weighted meansquare-error (wMSE): where x andx represent the noise-free and reconstructed signals, respectively. The ability to estimate macroparameters was assessed using the relative bias % computed as follows:  All animal experiments were conducted in compliance with French regulations and approved by the local animal ethics committee (APAFIS#6145).
The in vivo studies were carried out in female balb/cJ RJ mice (18.9 ± 1.3 g, from 8 to 10 weeks old) Half of the mice were xenograft subcutaneously on the right flank with 2.105 4T1 cells suspended in a solution of 100 μl of NaCl 0.9%. Imaging was performed 14 d after transplantation. Dynamic data were acquired with the Inveon TM PET/CT scanner (Siemens Healthcare, Erlangen, Germany). Mice were anesthetized for the duration of the imaging sequence by inhalation of isoflurane 2% and placed on a heated bed. After a 10 min transmission scan, the mice were injected into the caudal vein via a 28 ga catheter ( Image reconstruction was performed using CASToR software (Merlin et al 2018) with OSEM algorithm at 10 iterations and 24 subsets (Point-Spread Function (PSF) and post-smoothing). The reconstructed matrix was 200 × 200 × 170 voxels, with an in-plane voxel size of 0.398 × 0.398 mm 2 and a slice thickness of 0.796 mm. PET frames were reconstructed with increasing duration after injection: 6 × 5 s, 3 × 10 s, 4 × 60 s, and 300 s until the end of the scan or the next injection. PET images were reconstructed with time frames 4 × 15 s, 4 × 30 s, 3 × 60 s, 2 × 120 s, 5 × 240 s, 2 × 450 s, 3 × 900 s for a total duration of 5400 s. TACs were extracted from tumors and contralateral healthy tissues for cancerous mice, and from healthy tissues on both sides for healthy mice. IDIF was extracted from a region drawn within the heart ventricles.  Figure 4 shows the bias of Patlak K i for single-tracer compared to dual-tracer acquisitions, with various activity ratios. Bias was below 5% in single-tracer experiments and went from (1.1 ± 11.82)% to (46.15 ± 55.79)% in dual-tracer experiments as a function of activity ratio. An activity ratio 1/2 (using [ 18 F]FDG two times more than [ 18 F]FCH) or 1/1 (same quantity of radiotracers) provided the best compromise between bias and variance as compared to the single-tracer scenario.  Figure 5 shows Patlak K i values computed in healthy tissues and tumors for the single and dual-tracer acquisitions for [ 18 F]FCH and [ 18 F]FDG tracers. The kinetic parameters were computed with the same time windows as the ones used for the 90 min simulations. A significant variability for tumors K i estimations between mice (for both single and dual-tracer cases) was observed and K i was slightly overestimated in tumors for [ 18 F]FDG. However, they could still be well differentiated from K i in healthy tissues.

60 min scenario ([ 18 F]FDG+[ 18 F]FCH)
3.2.1. Single-tracer K i Figure 6 provides the bias of K i estimation for the single-tracer acquisition of 60 and 30 min. [ 18 F]FDG K i bias was around −4% for the 60 min scan and increased to 17% for the 30 min scans. Estimated [ 18 F]FCH K i values were less affected by the duration of the scan: it was maintained around 2% ± 7% for scans of 60 to 30 min. It was apparent that the K i estimation of [ 18 F]FCH was more robust in terms of duration, due to the fact that the kinetic of [ 18 F]FCH can be studied in the first 30 min of acquisition (Sutinen et al 2004, Schaefferkoetter et al 2017. Therefore, in the 60 min scenario, it seems more reliable to use [ 18 F]FCH as the second tracer. Seemingly, it is preferable to inject [ 18 F]FDG first since it requires 60 min to achieve favorable results.  For single-tracer, we observed that more reliable kinetic parameters were obtained when using TACs estimated from the SA coefficients compared to directly using noisy TACs (see appendix A.2). Figure 7 shows the wMSE values of IDIF fits for dual-tracer simulations with varying activity ratios between tracers, taking the single-tracer acquisitions as reference. The first tracer was set as [ 18 F]FDG, and the second as [ 18 F]FCH 30 min later with an additional scan time of 30 min for a total scan of 60 min. The wMSE was less than 8% for any case. For dual-tracer injections, the best result was for a 1/1 activity ratio with wMSE less than 2.5%. Figure 8 provides the bias of K i estimations for single-tracer simulations and for dual-tracer measurements with varying activity ratios. The single-tracer boxplots (60 min [ 18 F]FDG and 30 min [ 18 F]FCH) are displayed as references to evaluate the dual-tracer results. It was anticipated that larger quantities of tracer provide accuracy in K i estimation. The [ 18 F]FCH bias increased from (0.01 ± 13.85)% to (34.4 ± 107.12)%, whereas [ 18 F]FDG K i bias varied from (2.85 ± 33.99)% to (9.04 ± 19.76)% when the activity ratio ranged from 1/4 to 4/1. The best compromise between bias and variance was obtained for the 1/2 activity ratio.  For single and dual-tracer measurements, we also observed that applying Patlak to estimated TACs (from the SA coefficients) provided more robust results than directly using SA coefficients to compute the kinetic parameters (see appendix A.4).

Discussion
Dual-tracer acquisitions allow to eliminate physiological changes and external or internal movements that would occur between two independent scans separated by several days. The main goal of this study was to develop a generic approach to study the kinetic parameters of a single dual-tracer PET imaging acquisition with the least number of assumptions and parameters. Incorporating the general SA, arguments from Slawski and Hein (2013) and the NNLS algorithm, allowed to promote sparsity of the spectral coefficients with no requirement of any tuning parameter as opposed to the method proposed by Bell and colleagues (Bell et al 2017). The parameterization of the SA simply requires that the basis functions cover an appropriate range for both smooth and transient events of the tissue tracer concentration. The fitted spectral coefficients were finally used to recover the separated signals of each tracer to be able to perform more specific kinetic analysis on these separated signals. Moreover, SA did not only separate the tissue TACs of dual-tracer measurements, but it also acted as a denoising method for both single and dual-tracer acquisitions (appendix A.2).
In order to investigate the efficiency of the proposed approach, [ 18 F]FCH and [ 18 F]FDG signals were simulated using parameters extracted from our preclinical experiments and then contaminated with Gaussian noise to obtain realistic data similar to the clinical ones. According to Sutinen et al (2004), Schaefferkoetter et al (2017), it is recommended to compute kinetic parameters of [ 18 F]FCH in the first 30 min after injection. In the case of [ 18 F]FDG, it is more common to obtain kinetic parameters from the late frames (at least 30 min after injection) (Hu et al 2020). Therefore, inspired by Kadrmas et al (2013), the efficiency of the proposed technique was examined on simulations of 90 min ([ 18 F]FCH+[ 18 F]FDG) acquisitions with a 30 min delay between injections. Preclinical experiments on healthy and cancerous mice were performed with the same 90 min ([ 18 F]FCH+[ 18 F]FDG) scenario for validation. The number of animals involved in the study was not significant for statistical analysis. However, for both healthy and tumor tissues, similar K i values were obtained for single and dual-tracer scans. The healthy and tumor tissues could still be well distinguished with dual-tracer acquisitions.
With the purpose of using the method in clinical routine, we finally aimed at reducing the acquisition duration to 60 min. To select the injection order for this scenario, the K i bias for both tracers was compared to each other for scan durations of 60 and 30 min. Under this constraint, [ 18 F]FDG was selected to be injected before [ 18 F]FCH since it is required to have a 60 min duration acquisition. With a delay of 30 min between injections and 1/2 activity ratio, the method demonstrated the ability to well separate the two tracer signals from the dual-tracer acquisition with reasonable bias on K i values. Additional preclinical experiments on healthy rats were used to assess the robustness of the method in another context and with other tracers (see appendix A.1). AIF signals were also extracted from arterial blood samples. Estimated K i and V d values of [ 18 F]FDG and [ 11 C]Acetate for single and dual-tracer measurements seem coherent with each other. A bias was observed on [ 18 F]FDG K i values obtained with IDIF as compared to AIF, likely due to the spill-over effect when extracting the IDIF from heart ventricles. However, the variability of the results with IDIF was smaller than with AIF and the estimations were more precise. This enforces the idea of using IDIF for a wider application of the proposed generic method.
The main limitation of this work is the fact that the method is applied after tomographic reconstruction of the dual-tracer images so that single-tracer images are not available. In order to reconstruct separate images of each tracer, the next step of our work is to include the proposed method directly within the reconstruction task using nested optimization (Wang and Qi 2010) such as in Chalampalakis et al (2021).

Conclusion
We proposed a generic and robust method for analyzing dual-tracer PET acquisitions without requiring any tuning parameters. This approach could be easily used for any kind of dual-tracers couples so that to, firstly, derive the optimal injection framework and, secondly, solved the separation problem of the two radiotracers within this chosen framework.
128 × 128 voxels, with an in-plane voxel size of 0.38 × 0.38 mm and a slice thickness of 0.79 mm. The nominal in-plane resolution with fluorine 18 was 1.4 mm full width at half maximum in the FOV center (Bao et al 2009). PET frames were reconstructed with increasing duration after injection: 6 × 5 s, 3 × 10 s, 4 × 60 s, and 300 s until the end of the scan or the next injection. TACs were extracted from the left and right cortex using an atlas. IDIF were extracted from regions drawn in the heart ventricles.
In figure 9,  A.2. Impact of tissue TAC estimation using SA Figure 10 displays the bias of Patlak K i estimations for single-tracer acquisitions in two cases: using SA to estimate the tissue TAC or without it (i.e. using directly the noisy TAC). Comparing the results, it can be seen that the results with SA estimation are less biased and the variance is smaller. Indeed, SA estimation not only helps to separate the tissue TACS for dual-tracer acquisitions but is also able to denoise the tissue TACS in both cases of single and dual-tracer. In this work, the single-tracer results taken as a reference were obtained using the SA to denoise the TACs before performing the kinetic analysis. Acetate V d of single and dual-tracer measurements (b) with AIF and (d) with IDIF for each rat (one color per individual). Two K i s and V d s are estimated for each rat: one for the right cortex which is plotted by a plus (+) and the one for the left cortex is displayed by a cross (×).
A.3. Impact of noise on K i estimation In our simulations, we fixed the intensity of noise using ν = 0.35 to empirically mimic realistic IDIF and TACs. Figure 11 shows the bias of K i estimations using the Patlak plot for different intensities of noise (ν = {0, 0.35, 0.7} in the case of dual-tracer 60 min acquisitions with activity ratio 1/2. It can be seen that the bias distributions of K i estimations for both radiotracers become larger with increasing noise levels while retaining meaningful averaged bias values. This figure confirms the robustness of the proposed strategy with respect to noise intensity.
A.4. Bias of K i using SA The net uptake of the tracer in the tissues K i in an irreversible process using SA is given as: figure 12 displays the bias of K i estimations using SA model (equation (8)). Compared to the results of K i from Patlak plot (figure 8), the latter are more robust and the error is smaller in both single and dual-tracer settings. As it can be seen in figure 12, it is not possible to estimate K i using SA in the case of 30 min single-tracer [ 18 F]FCH, whereas Patlak plot demonstrates accurate and precise K i estimations in this situation ( figure 8). Therefore, it is more reasonable to estimate K i using the Patlak plot where the SA is used only to separate and denoise the tissue TACs.
A.5. Simulation of the tissue TAC using the 2TCM In 2TCM (illustrated in figure 13), compartments are connected with four rate constant K 1 , k 2 , k 3 and k 4 . For an irreversible process, we have k 4 = 0. The tissue TAC is the sum of the two tissue compartments as  C(t) = C 1 (t) + C 2 (t). The derivative of tissue compartments C 1 (t) and C 2 (t) with respect to t can be computed as: Therefore, C 1 (t) can be written as: . 11 Using the trapezoidal rule, the integral of X can be written as Using this rule and (11), we can write: and C 1 (t − Δt) respectively represent the integral and value of C 1 over the previous frame. From (10), C 2 (t) can then be computed as: Eventually, in order to compute the tissue TAC in 2TMC, (11) and (10) are substituted in C(t) = C 1 (t) + C 2 (t).