Quantitation of Multiple Injection Dynamic PET Scans: An Investigation of the Benefits of Pooling Data from Separate Scans when Mapping Kinetics

Multiple injection dynamic positron emission tomography (PET) scanning is used in the clinical management of certain groups of patients and in medical research. The analysis of these studies can be approached in two ways: (i) separate analysis of data from individual tracer injections, or (ii), concatenate/pool data from separate injections and carry out a combined analysis. The simplicity of separate analysis has some practical appeal but may not be statistically efficient. We use a linear model framework associated with a kinetic mapping scheme to develop a simplified theoretical understanding of separate and combined analysis. The theoretical framework is explored numerically using both 1-D and 2-D simulation models. These studies are motivated by the breast cancer flow-metabolism mismatch studies involving 15O-Water (H2O) and 18F-Fluorodeoxyglucose (FDG) and repeat 15O-H2O injections used in brain activation investigations. Numerical results are found to be substantially in line with the simple theoretical analysis: mean square error (MSE) characteristics of alternative methods are well described by factors involving the local voxel-level resolution of the imaging data, the relative activities of the individual scans and the number of separate injections involved. While voxel-level resolution has dependence on scan dose, after adjustment for this effect, the impact of a combined analysis is understood in simple terms associated with the linear model used for kinetic mapping. This is true for both data reconstructed by direct filtered backprojection (FBP) or iterative maximum likelihood (ML). The proposed analysis has potential to be applied to the emerging long axial field-of-view PET scanners.


Introduction
Positron emission tomography (PET) scanning is an important diagnostic imaging technique used in the management of patients with cancer and other diseases, as well as in medical research. In clinical settings, PET is typically used to create a single 3-D static image of the distribution of the tracer atoms in the field of view in a time-frame of interest, after the tracer has been administered. But PET can also be used to obtain a temporal sequence of scans after radiotracer injection. Such dynamic PET scans give the possibility to analyze the transport and retention of tracer in order to recover more detailed metabolic information (Cunningham and Jones;1993;1993). Most dynamic PET studies involve the use of a single radiotracer, such as 15 O-Water (H 2 O) for blood flow 1983), 18 F-Fluorodeoxyglucose (FDG) for glucose metabolism 1979) or 11 Cverapamil for P-glycoprotein (P-gp) activity (Deo et al.;. Studies in which multiple radiotracers are used in the same subject in the same imaging session have the potential to provide more comprehensive profiles of tissue function , in-vivo. Examples include: repeat 15 O-H 2 O imaging for brain activation studies (Beason-Held et al.;, 15 O-H 2 O and 18 F-FDG for breast cancer , 13 N-NH3 and 18 F-FDG for myocardial viability (Yamagishi et al.;Dilsizian et al.;, 15 O-H 2 O and 11 C-Verapamil for assessment of multi-drug resistance (Eyal et al.;2009).
While data from multiple radiotracer scans could be assembled post-hoc from separate imaging sessions, acquiring the data in a single imaging session has several advantages. (i) It can reduce issues of anatomic co-registration. This is important for comparing the facets of biology related to multiple tracers, especially for targets that may move between studies, for example, breast (Fowler et al.; and ovarian (Makvandi et al.; cancers comparing FDG PET/CT staging studies to therapeutic biomarkers. No doubt longer separation of injections can lead to more patient motions. (ii) It can reduce potential variation in the physiologic state of the subject from one session to the next 1982). This is especially true for metabolic studies using H 2 O and FDG or Glutamine and FDG (Zhou et al.;. (iii) Tracer parameter estimation for multi-tracers may provide more reliable kinetic analysis due to common dependence of parameters: tracer delivery on flow or tracer clearance on plasma volume and hydration, etc.
1998; Rust and Kadrmas;2005;Zhang et al.;. However, whenever the temporal separation between different tracer injections is much longer than the radiotracer half-life or there is a simple mechanism to correct for spillover 1982;, 2001Black et al.;2009;El Fakhri et al.;, analysis can be approached in one of two ways: (i) studies can be processed separately to recover kinetic information corresponding to individual tracers ( Figure 1A), or (ii), studies can be concatenated temporally and a combined analysis ( Figure 1B) of the resulting data carried out. This paper develops a simplified theoretical framework for understanding the relative performance of separate and combined analysis of multi-tracer PET studies. The approach is based on a type of linear modelling framework for mapping kinetics 1993; -details in section 2. The theoretical analysis suggests that after controlling for voxel-level imaging accuracy (resolution), the relative performance of separate and combine kinetic mapping can be understood in simplified terms. A series of numerical experiments, motivated by flow-metabolism mismatch and repeat brain activation studies, are used to directly assess mean square error (MSE) characteristics of combined and separate analysis techniques. These studies consider both filtered backprojection (FBP) and maximum likelihood (ML) reconstructed data as well as a number of important factors including, overall dose, the relative dose of individual tracers and the number of injections involved. Data from numerical experiments are subjected to careful statistical analysis. This helps to clarify the main findings. In flow-metabolism studies, combined analysis improves MSE performance by between 8% and 30%, depending on the reconstruction method used. In the repeat injection studies, MSE improvements increase linearly with the number of injections involved. This is in line with the proposed theory.
The outline of the paper is as follows: The kinetic mapping approach and the associated analysis of the impact of data combination are developed in section 2. An illustrative practical example from a flow-metabolism mismatch study in a breast cancer patient is given in section 3. Section 4 describes a series of numerical studies based on 1-D and 2-D scanning models. Results are presented in section 5. Section 6 presents a detailed statistical analysis of simulation results. The paper concludes with discussion.

Theory
The focus is on dynamic PET studies in which it can be assumed that the tracer's interaction with tissue is substantially linear and time-invariant. This assumption is valid for a large class of tracers used in PET, particularly in a cancer setting. The interpretation of timecourse data from dynamic studies can be approached using classical indicator dilution theory (Meier and Zierler;1954). This theory expresses the measured time-course as a convolution between the arterial input function (AIF) and the so-called tissue residue or impulse response function. The tissue residue is a necessarily monotone decreasing function. Widely used compartmental models (Kety and Schmidt;1948;Phelps et al.;1979), typically approximate the residue as a positive linear combination of 1 or 2 mono-exponential functions. The spectral method (Cunningham and Jones;1993) generalizes this to allow the positive combination of potentially many mono-exponentials. In statistical terms, the residue may be viewed as a life-table for tracer atoms introduced to the tissue at time zero via the arterial supply -indeed this is a key part of (Meier and Zierler;1954). As life-tables can be represented in terms of the underlying distribution (density/histogram) of travel times associated with the physiologic/metabolic journey of individual tracer atoms in the tissue, any residue modelling assumption carries with it a specification for the distribution of such travel times within the volume of tissue being analyzed 2009). Compartmental and spectral approximation requires that travel-time densities be spiked at time zero and be monotone decreasing. However, strict monotonicity of travel-time densities may not always be reasonable 1997). One might imagine this may only be practically important in studies with high-frequency temporal sampling -perhaps only in the context of the emerging PET scanning technologies (Badawi et al.;Karp et al.;2020). However 2009) used circa-1996 PET scanner data involving relatively crude temporal sampling to demonstrate significant statistical deficiencies in the ability of the standard 2-compartmental model to correctly represent the regional cerebral time-course of PET-FDG data in highly homogeneous brain regions of normal subjects. The anomaly was identified by considering a more flexible non-parametric analysis of the residue function, with suitable cross-validation adjustment for potential overfitting relative to the 2compartmental model. In light of this finding, we have been motivated to employ a nonparametric representation of the tissue residue, particularly when analyzing PET time-course data corresponding to large and potentially heterogeneous regions of interest. The concept has been incorporated into a non-parametric residue mapping (NPRM) methodology that we have found useful for parametric imaging of dynamic data from PET imaging of cancer 1993;. Figure 2 shows a schematic. Before we consider the topic of combined analysis of multiple-injection studies, we elaborate on important details of the NPRM approach. Theoretical points are illustrated numerically in the context of kinetic mapping PET-FDG data.

Voxel-level Kinetic Mapping by the NPRM Approach
There is extensive literature on techniques for mapping kinetics from dynamic PET datasee 2020) for a recent review including many references. A broad family of techniques, including (Cunningham and Jones;1993;1993), involves voxellevel representation of the tissue residue as a (positive) linear combination of suitable basis functions. Mono-exponential basis residues are used in (Cunningham and Jones;1993); more elaborate basis elements, incorporating potentially complex blood-tissue exchange models are proposed in 1993). NPRM technique is a variation on the latter approach in which basis elements (known as sub-TACs) are represented in terms of flexible model-free residues. A full account of the method with several cancer imaging applications is provided in . In the current implementation of NPRM, sub-TACs are constructed by a 2-step data-adaptive process. We elaborate on these steps and indicate how derived kinetic are mapped.
Step 1 applies a standard recursive hierarchical clustering scheme to partition the image volume into a set of S segments or clusters, with the property that the time course data for voxels in the same segment are each approximately proportional to the mean time-course for that segment. Typically the number of segments needed for this is on the order of 100-200. The mean PET-measured time-course data for segments (z s ) together with their standard deviations (σ s ) are used to produce a reduced data set D S = z sl , σ sl , s = 1, 2, …, S, l = 1, 2, …, N T . Here N T is the number of time-frame of PET data acquisition. Next, the mean time-course data for each segment is modelled to create an associated set of modelled elements ℳ = μ s , s = 1, 2…, S . In the case of a single-injection, a generic time-course (z) in D S is modelled (μ) by one of: (a) a scaled (venous) injectionsite signal (C IV ), (b) a non-parametric residue function (R), or (c), a scaled non-parametric distribution function (F) representing the total number of tracer atoms no longer in bloodtissue exchange as a function of time. This last case is to accommodate data in which the bladder is in the field of view. Thus μ l = A ⋅ C IV t l − Δ ; μ l = R ⋆ C p t l − Δ ; μ l = F ⋆ C p t l − Δ (1) for l = 1, 2, …, N T . Here C p is the time-course of the tracer in the arterial blood (AIF). Δ is an optimized delay factor. R and F are non-parametrically specified as piece-wise linear functions and estimated using constrained weighted least squares with weights 1/σ l 2 . It is important to appreciate that, apart from Δ, all unknowns in (1) are linear, so the computation of the optimal fit, for any fixed Δ, is evaluated by quadratic programming (QP). A grid search is used to determine the optimal delay. Choices other than the residue case are only selected if they significantly improve on the residue fit. Modelled elements are supplemented to include time-courses corresponding to the AIF and a further time-course corresponding to the integral of the AIF -the AIF exactly models a spike residue; the integral of the AIF corresponds to a constant residue. In deference to (Patlak and Blasberg;1985), the latter component is referred to as the Patlak-element. With multiple injections an elaboration of (1) is used -details of this are indicated below.
Step 2 of the basis selection procedure takes the full collection of modelled time-course patterns and applies a cross-validation guided backwards elimination procedure to identify a subset of these with the property that positive linear combinations of these patterns adequately fit the reduced dataset D S . This leads to a final set of basis elements given by a N T × K matrix X = μ 1 , …, μ K . These basis elements (μ) are referred to as sub-TACs and K is their total number. Given the manner in which the segmentation process is implemented, voxel-level time-course data (z i ) can be expected to be well-approximated by a positive linear combination of sub-TACs. NPRM fits voxel-level data by the linear model and approximates the voxel-level residue using the components that correspond to a residue process (1) where R t ⋅ , x i is the residue for the i'th voxel. Here K r is the indices of the subset of sub-TACs that are specified by a residue in (1). QP is used to evaluate α i .

Residue Decomposition and Mapping of Kinetics:
Suppose a generic residue function R can be measured over a time interval [0, T E ] -the PET scanning time-frame. Let T B for 0 < T B < T E represent a realistic upper-bound for the large-vessel traveltime of tracer atoms -a T B value of 15 seconds seems physiologically reasonable, for most PET tracers, including FDG and H 2 O. Following , the residue can be written as a sum of vascular, in-distribution and extracted elements These assess, large vessel and small vessel distribution volumes (V B , V D ), non-large vessel flow (K D ) and retention or apparent flux (K i ) at the end of scanning period. Note, that all of these quantities are linear functions of the residue. Tracer extraction is measured as K i /K 1 where K 1 = K D + K i represents the overall flow. By the central volume theorem (Meier and Zierler;1954), mean transit time (MTT) of non-extracted tracer is measured by a ratio of volume to flow. The measure may be adjusted to account for an overall delay (Δ) in the arrival of tracer atoms to the tissue region. This gives MTT = V D /K D + Δ.
In the NPRM procedure voxel-level kinetic parameters are based on the estimated residue in (2). This leads to mapped values for (V B , V D , K D , where Δ is the flow-weighted average of the delays (Δ k ) associated with residue components at the voxel being mapped, i.e. at voxel i, Δ = ∑ k ∈ K r α ik K Dk Δ k ∑ k ∈ K r α ik K Dk with K Dk the flow for R k .
Note (V B , V D , K D , K i ) are linear and (MTT, K i /K 1 ) are (smooth) non-linear functions of αcoefficients in (2).

Why is the NPRM Approach Reasonable?-
The NPRM approach has three important features: (i) In diverse tissue environments the formal theoretical support for compartmental models may be limited (even from an in-vitro stand-point), especially with more novel PET tracers. In this context a non-parametric approach has improved flexibility to adapt to the true (unknown) kinetics. (ii) NPRM can accommodate significant variations in the arrival of tracer to different parts of the field of view. This is accomplished by the inclusion of delays in the modelling of sub-TACs. (iii) NPRM gives an ability to adjust for potential artifacts associated with the bladder or the injection site within the scanner field of view. This is not important in the data presented in this paper but has relevance in abdominal scanning and in cases where there may be enhanced temporal scanning resolution.
An alternative non-parametric approach, based on a truncated singular value decomposition (SVD), has been used to map perfusion parameters with MR and CT data (Østergaard et al.;1996;Konstas et al.;2009). Unlike the method used here, the SVD approach does not constrain the target residue to be either positive or monotone decreasing. Also, it does not have a mechanism to address heterogeneity voxel-level data or indeed to accommodate artefacts associated with the retention of tracer in the bladder or signals associated with a venous injection site in the field of view. Further examination would be helpful in getting a detailed clarification of the performance of the truncated SVD approach relative to the NPRM technique. NPRM mapping would also seem to have potential for application to perfusion imaging with MR and CT. In addition, the approach may also have a role in quantitative analysis of digital subtraction angiography imaging (DSA) data -see (Liang et al.; A small 1-D PET-FDG numerical simulation study comparing NPRM with classic voxel-byvoxel compartmental analysis (Kety and Schmidt;1948;Phelps et al.;1979) is implemented.
Details of the scanning model and source distribution are as specified in section 4. The source distribution is specified as a linear combination of specified sub-TACs (μ k ) -like (2). In one case (I), source sub-TACs correspond to non-parametric residues derived from regional analysis of breast cancer imaging data; in the second case (II) source sub-TACs are replaced by the best fitting one-compartment (1C) model (Kety and Schmidt;1948) to the real data. The table reports the relative root mean square error (RMSE) performance of mapping the six summary kinetic parameters described above, when voxel-level residues are estimated using either the NPRM technique or by application of a 2-compartment (2C) model 1979) that includes adjustment for the fractional blood volume and tracer delay 2009). The model function (η) is: where C p is the AIF and f b accounts for the contribution from the vasculature. The unknown parameters are β= (Δ, f b , A 1 , A 2 , λ 1 , λ 2 ) -Δ is unconstrained but the other components of β are constrained to be positive.   Table 1 presents the comparisons between the NPRM and 2C methods when the source sub-TACs are generated by a compartmental form (case II) and also when the sub-TACs follow a non-parametric form (case I). Results show substantially enhanced performance of NPRM when the source is not compartmental and a more similar performance when the source is compartmental. Remarkably however, perhaps due to spatial mixing associated with data reconstruction, the 2C method does not perform as well for the V B and V D parameters.
Indeed NPRM is seen to be strongly competitive with the 2C mapping scheme when the underlying source is specified in terms of a compartment model structure. At its worst, in the case of K i /K 1 , the relative RMSE efficiency of NPRM mapping is 84%. This result is quite typical.

Local Error Characteristics of NPRM Mapping:
Although the voxel-level model in (2) used by NPRM has a linear form, the sub-TACs involved, X = μ 1 , …, μ K , are derived by a complex data-adaptive process. Hence the overall error in voxel-level kinetics is a function of the uncertainties in the estimation of α-coefficients and the uncertainties in the construction of sub-TACs. Intuitively, because the sub-TACs are constructed by modelling average time-course data for segments, one might expect that the uncertainties associated with the sub-TACs have a minor impact on the error on NPRM voxel-level kinetics. In order to justify this intuition, we report on another simulation. Here we used a standard 2D PET-FDG simulation in which voxel-level kinetics were assessed using both known and unknown data-estimated sub-TACs. The 2-D context is used as the segmentation process with 2-D data gives rise to more realistic subsets of data defining each sub-TAC. The scanning model and source structure are as described in section 4. The statistical relation between the RMSE (averaged across voxels) using known and unknown sub-TACs are shown in Table 2. This study only considered a single (realistic) dose and was restricted to a modest number of replications (N R = 50) -standard errors of simulation estimated quantities of interest do not suggest we need any more. The results show that kinetic mapping errors are substantially the same regardless of whether or not sub-TACs are known. Thus when analyzing voxel-level NPRM kinetic mapping errors, it seems sufficient to focus on error associated with determination of voxel-level α-coefficients conditional on known model sub-TACs. We will adopt this approach in our analysis of the impact of data combination.
We now focus on what controls NPRM kinetic mapping error. A simplistic linear model analysis (Seber;2009) of the NPRM model (2) would suggest that the RMSE of the kinetic summaries should substantially scale with the level of uncertainty in the reconstructed data. Thus for a given voxel (i) the error in an estimate of kinetic variable, θ ip , denoted RMSE ip NPRM should have a strong relation to the corresponding root-mean-square deviation between the reconstructed time-course data (z il ) and the true voxel-level time course We examined this numerically using the simulation data used to produce Table 1. The statistical analysis, which incorporated data from a range of dose settings, finds log RMSE ip /θ ip ≈ α p + β p log σ i /θ ip (7) well-describes the relation between voxel-level errors in kinetics and the local reconstruction accuracy of the input data (σ i ). Model R 2 -values and β-coefficients are given in Table 3.
Note that (7) suggests formal mean square error convergence of the NPRM procedure -as the reconstruction error diminishes so too does the error in determination of kinetics. The rate of convergence, measured by the β p coefficient, is seen to be close to 1 and suggests that standard parametric inference asymptotics apply to this situation (Van der Vaart; 2000). A similarly strong relation is found between the modelling error associated with the linear In simulations this is important as it allows us to use modelling error as a surrogate metric for the analysis of errors in individual kinetic variables.

Multiple Injections and their Combined or Separated Analysis by NPRM
With multiple injections the model in (1) is elaborated to account for the effect of individual injections and the separate delay factors that may be associated with measurement of the tracer AIFs. This gives Here the index j relates to individual injections. As with (1), the residue and outflow functions can either be modelled using suitable parametric forms 1993;Spence et al.; or with non-parametric formulations as used by (Lan et al.;. Note the non-parametric approach leads to an optimization problem, in which only the delay factors Δ j enter in a non-linear fashion. This simplifies computation. Our interest is in the case where it is possible to consider the separate analysis of kinetics of individual injections. In this setting, because there is negligible temporal overlap (spillover) between the separate injections, modelling of sub-TACs in the NPRM voxel model can be accomplished by separate application of the simple form in (1). Regardless of whether there is overlap or not, voxel-level kinetic parameters, say {θ p , p = 1,2, … P}, are expressed as simple functions of the α-coefficients in (2). Thus θ ip = g p α i where g p is a (smooth) function defined using the estimated residue functions and delay factors for separate injections.
If data are combined, coefficients are obtained based on the full time-course for the measured data and the sub-TACs. Separate analysis recovers kinetics based on the shorter sub-TAC time-course corresponding to individual injections. Now suppose we are interested in a kinetic summary variable corresponding to the j'th tracer injection -say θ = g(α). The separate analysis estimate, θ j = g α j , is based on fitting (2) using only the j-injection timecourse; the combined estimate, θ C = g α C , is based on fitting the full time-course data. Thus at the i'th voxel we would map the kinetic parameter associated with g using one of Intuitively, we expect that the combined data estimator will be preferred but it is helpful to have a quantitative appreciation of the factors that might be important. The results reported in the paper address this issue via numerical simulation. Before that, we present a simplified theoretical analysis that may give some intuition. Table 3 suggests that regular parametric statistical inference theory associated with Gaussian approximation applies to the NPRM model coefficient estimates and the derived kinetic variables. In light of this it seems reasonable to use Gaussian approximation to develop insight into the potential benefits of combined versus separate analysis for multi-injection studies. By Gaussian approximation we get the following.

An Approximate Theoretical Analysis of Combined versus Separated Analysis
Result 2.1.-If the conditions for Gaussian approximation hold, the combined analysis estimator always has lower variance than the separate injection estimator. For a given kinetic parameter (θ), the ratio of the MSE of the combined analysis estimate, based on J injections, to the MSE of the separate analysis estimate, using only study j, is given by the formula where σ j 2 is the variance of the kinetic parameter based on data from injection j.
The Appendix provides a theoretical justification for this result. In replicate studies, individual scan variances (σ j ) are identical so the MSE ratio is 1/J, i.e. the percent improvement in MSE with the combined analysis is proportional to (J − 1). While the combined analysis will always produce better estimates, the amount of improvement will depend on the reliability of the other scans involved. (8) indicates that the reliability of kinetics for individual scans kinetics is a function of the associated source model error. We would anticipate that the study dose and relative sharpness of the AIF, which impact voxellevel model error, would play a role. Numerical studies will explore this in some detail.

Illustration with a Flow-Metabolism Study in a Breast Cancer Patient
The data come from a study conducted on a breast cancer patient prior to surgical resection and prior to scheduled neoadjuvant chemotherapy. PET scanning involved dynamic imaging with 15 O-H 2 O and 18 F-FDG in the same session. 15 O-H 2 O is a tracer used to measure blood flow and perfusion -the gold-standard technique for in-vivo measurement of such information. 18 F-FDG is the most widely used clinical PET tracer. FDG measures information about tissue glucose utilization. It is a key marker used in the diagnosis and clinical management of several cancers. Details of PET dynamic imaging protocols for H 2 O-FDG scanning in breast cancer and reports on the prognostic utility of the derived information are given in the earlier reports -see , 2003. For the data we present here, 1302 MBq of 15 OH 2 O in a 1-4 mL volume was injected as a bolus and 318 MBq of 18 F-FDG in 7-10 mL volume was injected over 2 minutes with a constant infusion pump. Data were acquired on a GE-Advance scanner using a plane-by-plane FBP reconstruction algorithm with corrections for attenuation, scatter, deadtime and random events. The 4-D PET data set consists of an imaging volume with N =  Figure 3. Coronal planes of metabolic parameters demonstrate increased flow and metabolism in the tumour. The flow-metabolism mismatch is also seen to be highly elevated in the tumor region -an elevated mismatch has been shown to be associated with poor response to chemo-therapy and early cancer relapse . In qualitative terms the flow and mismatch maps recovered by separate analysis appear more noisy, the flux maps are not much different. To explore these differences more carefully, two regions-of-interests (ROIs)normal breast (335 voxels) and tumour (322 voxels) are extracted, as shown in Figure 3. The distribution of metabolic parameters in these two ROIs is presented in boxplots in Figure 3.
We can see the flow distributions in tumour and normal breast are more variable for the separate analysis, the standard deviations (SD) of flow from combined and separate analysis in 2 ROIs are 0.04 and 0.07 (tumour), 0.014 and 0.018 (normal breast) respectively. The additional variability is on the order of 75% and 28.6% in tumour and normal breast. For flux, the standard deviations with the two approaches are 0.006 and 0.005 (tumour), 0.0004 and 0.0003 (normal breast). These differences are small. The variability of mismatch from separate analysis is much bigger than the combined analysis in the tumour region -the situation is somewhat reversed in the normal region but this could also be due to an underestimation of the normal tissue flow in the separate analysis. In light of these results, we are motivated to develop a better understanding of the relative statistical behaviour of combined and separate analysis of this and similarly structured multiple injection studies.

Numerical Experiments
We conducted two sets of experiments, one related to the H 2 O-FDG study in section 3 and a second one relating to repeat H 2 O activation studies. A schematic for the structure of the simulation process is provided in Figure 4. Below we provide details of the setup of the 1-D and 2-D simulation models used, the specification for the source distributions in the two examples and, finally, the procedure used to compare the performance of the combined and separate analysis techniques. FBP and ML reconstructions are both considered in 1D simulation studies, but for computational efficiency we restrict the H 2 O-FDG setup to only FBP reconstruction in the 2-D case.

Scanning Models
In the 2-D case a standard PET scanning model involving Poisson sampling of a discretized (attenuated) parallel-beam Radon transform of the source is used (Kak et al.;Natterer;2001). The imaging domain is the unit square, discretized to an array of dimension 128 × 128, and the projection domain is the region [ − 2, 2] × [0, π], discretized to a 183 × 181 sinogram array of distances and angles. For computational reasons numerical experiments in 2-D only consider direct FBP reconstruction. In 1-D a less familiar Poisson deconvolution model is employed. The model is used because the computational efficiencies enable us to carry out a sufficient number of replicate studies using both direct and iteratively reconstructed data. Full details of the model are given in (O'Sullivan and Roy Choudhury; 2001). The observed data over a given time-frame is a realization from a projected and scaled source. The projection process is given by the product of a fixed attenuation profile (a > 0) and the discrete convolution between the scaled source and a kernel κ β (β > 0). The kernel acts to smooth the source -its discrete Fourier transform is given by κ ν β = ν −β for ν = ± 1, ± 2…N /2.
If y is a realization from the projected and scaled source then the direct reconstruction of y is simply obtained where (y/a) is the discrete Fourier transform of the attenuation corrected count data, ℱ −1 represents the inverse transform and τ is the expected total counts (dose) used for scaling.
Smoothing is achieved by convolving the raw reconstruction z* with a discrete Gaussian kernel. The smoothed reconstruction is an analogue of an FBP reconstruction used in PET. ML reconstructions are readily obtained in this setting. Standard iteratively re-weighted least squares (IRLS) (with positivity constraints) apply. The convolution structure allows such iterates to be easily evaluated. Similar to FBP, the raw ML estimator can be smoothed to obtain an analogue of the expectation-maximization (EM) reconstruction methods used in PET. The simple convolution model captures the essential estimation complexity of PET. This is because the reconstruction filter acts like a fractional derivative. By choice of β, it is possible to adapt the model so that the behaviour of the reconstruction MSE as a function of count rate, is matched to 2-D PET. Experimentation by 2021) with the breast FDG data in section 3 found that a value of β = 1:35 in the 1-D scanning model worked well. This is the setting we use in the studies reported here.

Source Patterns
The true (target) source is expressed as a linear model with K components.
λ(x, t) = α 1 (x)μ 1 (t) + α 2 (x)μ 2 (t) + … + α K (x)μ K (t) Gu et al. Page 12 In the 2-D case we only consider the H 2 O-FDG studies. The configuration is matched to the real data presented in section 3 -see Figure 5. Spatial patterns correspond to a central slice through the tumor region; the temporal patterns including the AIF come from the full data set.
The source specification is more abstract in the 1-D case. Figure 6 shows the spatial pattern of the α-coefficients and how they are transformed into the measurement domain. This aspect is the same for H 2 O-FDG and the repeat H 2 O studies. The shape of the AIFs for H 2 O and FDG match those in the breast cancer data. The time-course patterns for the H 2 O-FDG are also adapted from the breast cancer analysis; in the repeat H 2 O studies, a 1compartmental Kety-Schmidt (Kety and Schmidt;1948) model is used to create patterns relating to regions in a normal brain -parameters are given in Figure 6. Figure 7

Evaluation of Performance
Our studies use the model MSE characteristics of the estimated source, as a surrogate for understanding the MSE characteristics of estimates of the mapped kinetics. The motivation for this comes from Table 3 which showed a setting where MSE characteristics of mapped kinetics are proportional to the MSE reliability of the underlying source estimate. As indicated in Figure 1 the overall scan period, T, is made up of a set of J non-overlapping time intervals corresponding to the scanning periods for the separate tracer injections under consideration, T = T 1 ∪T 2 ⋯∪T J . If α C is the combined estimate coefficient, the corresponding estimated source is λ C (x, t) = ∑ k = 1 K α k C (x)μ k (t) and the squared error assessment is where D j and λ j are the durations and maximum intensities for the j′th scan. Integrals are evaluated by simple discretization. For the separate scan case, the estimate of the source for the j'th time-interval is based on the corresponding coefficients (α j ) -i.e. λ S (x, t) = ∑ k = 1 K α k j (x)μ k (t) for t ∈ T j . In this way an overall assessment for the separate scan analysis is obtained as Se λ S , λ . With simulated data, we evaluate the error of combined and separated scan estimators for study replicates and by averaging an MSE assessment for the source error is obtained. This gives MSE C = E Se λ C , λ and MSE S = E Se λ S , λ Relative percent improvement is assessed by

Results
We begin by presenting some sample results for the H 2 O-FDG and repeat H 2 O studies. This is followed by detailed MSE comparisons, including some associated plots.

Sample Metabolic Parameters Estimates
In the breast cancer study, H 2 O and FDG are used to measure blood flow and glucose metabolism, respectively. Metabolic parameter comparisons for a single replicate at the middle dose and dose-ratio are presented in Figure 9. True parameters are represented by black lines. Estimated parameters from combined and separate analysis are shown using dashed red and blue lines separately. We can see metabolic estimates from combined analysis are closer to the true parameters and show more advantages in H 2 O study. Metabolism estimated from these two approaches are both close to the true values in FDG study and difference can be negligible. The quantitation of H 2 O can get more benefits from the combined analysis and these results are also consistent with findings from the real dataset in section 3. Similarly, we also show one sample metabolic parameters -flow and blood volume (V B ) at the middle dose in a repeat H 2 O study (H 2 O*) in Figure 9. The combined analysis is seen to achieve improvements in both situations.

H 2 O-FDG Studies:
The MSE in (14) as the function of the dose and dose ratio from combined and separated approach in 1-D and 2-D simulations are presented in Table 4. Figure 10 presents the information graphically. MSEs from combined and separated analysis decrease with the dose and increase with relative dose ratio. In 1-D simulation studies, the improvements are 20.85% and 11.65% from FBP and ML reconstruction at the reference dose and middle dose-ratio level -matched to the real dataset and improvement is smaller based on the ML reconstructed data. Benefits by combined approach in 2-D simulation with FBP reconstruction is similar to 1D FBP performance.

Repeat H 2 O Studies-
The quantitative improvements for different numbers of injections at seven dose levels from two reconstruction algorithms are presented in Table 5.
We can see, the number of injections has a strong effect. See also Figure 11. These results show lower doses and more injections have higher MSEs. Dramatically smaller improvements are seen with ML. This may be because the bandwidth of smoothing process used in reconstruction is adapted based on the total update scan so there is a more limited impact on individual time frames; the ML procedure constrains the reconstruction to be positive and in doing so implicity smooths the voxel-level time-course data.

Statistical Analysis of Simulation Results
The highly structured nature of the MSE patterns in Figure 10 and 11, invites a more refined interpretation.

MSE values scaled by Image Data Accuracy
The overall accuracy of the input data is assessed by σ z (7). While we can scale MSE values by this, a simpler (and statistically equivalent) procedure is to scale by the voxel-by-voxel error in the total-uptake reconstruction. We label this σ z • 2 . For the H 2 O-FDG study, the scaled values are in Figure 12.
These data are well described by the model log MSE/σ z • 2 a 0 + a 1 I + a 2 DR (16) where I is the method indicator for combined or separate analysis and DR is the dose ratio. Table 6 reports estimates values for these coefficients. Note a 1 corresponds to about 21% (e 0.19 − 1) and 12% (e 0.11 − 1) improvements at middle dose-ratio with the FBP and ML reconstruction in 1-D simulation studies. The results from 2D simulation have the similar performance.
A similar analysis of the scaled MSE data from the repeat H 2 O study, leads to a model with an interaction term: log MSE/σ z • 2 β 0 + β 1 I + β 2 No . Studies + β 3 I × No . Studies The scaled MSE values are in Figure 13 with the model results in Table 7. The improvements of combined analysis in the repeat H 2 O studies with two injections are 55% and 15% for the FBP and ML reconstructed data, respectively. The impact of dose, although significant, is not substantial in comparison to the overall impact of dose on MSE.

Analysis of the Uptake MSE as a function of Dose
Analysis of MSE of total uptake from FBP and ML reconstruction as a function of dose gives a 2-factor statistical model: where τ is dose level, M is dose ratio (DR) for H 2 O-FDG study and number of injections for repeat H 2 O studies, respectively. The estimated coefficients γ 0 , γ τ and γ M are shown in Table 8. We can see γ τ is similar in different studies and reconstruction methods. The Gu et al. Page 15 Phys Med Biol. Author manuscript; available in PMC 2021 July 16.

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript intercept -γ 0 in FBP is always bigger than ML, which means uptake based on ML has smaller errors. Table 8 shows that the main effect of dose is described by τ −.42 -this is consistent with previous asymptotic analyses of ML and FBP in PET, e.g. (O'Sullivan;1995). The fact that the rate of convergence is so similar between the 1-D and 2-D simulation models is a good indicator that the 1-D PET framework captures the estimation complexity of 2-D PET.
The rate of convergence analysis allows us to develop a relation between percent improvements obtained by combined data analysis and the equivalent precent increase in dose that would be needed to ensure that the separate data matched the benefits of combined analysis. Figure 14 presents this equivalence. In light of Table 8, which gives a dose rate coefficient on the order of 0.42 for ML and FBP, a 10% improvement in MSE obtained by combined data analysis would require a 25% increase in dose if separate data analysis methods were used.

Discussion
We have presented a combination of theoretical intuition, real data and numerical simulation to examine the benefits of combined NPRM analysis of multiple injection studies in PET. This work has been fully motivated by practical imaging protocols. The results demonstrate that there are clear gains in MSE performance by a combined analysis approach. With H 2 O-FDG studies the MSE improvement is around 21% on the same order of the regional variance comparisons found in section 3. For the repeat H 2 O studies, improvements increase in accordance with the factor (J − 1) indicated by the theory in section 2. In general improvements are seen to be a function of intra-injection dose-ratio but are substantially independent of dose for FBP reconstructed data. In the ML case, the impact of dose is apparent with combined analysis leading to greater improvements at higher doses. Even though the theoretical analysis looks at the asymptotic situation in which variance is the main contributor to the error, the predictions from that theory provide good guidance on the nature of benefits that would be associated with combined analysis.
Our work has exclusively focused on parameter mean square error performance (MSE and RMSE). In kinetic mapping, the average accuracy of estimation of the target parameters, measured by the RMSE, is a natural metric and is used extensively in quantifying estimation error in many settings. Nevertheless, a number of other metrics could also be of interest. Data-fit or predictive error, with suitable adjustment for model flexibility, might also be considered. In addition, criteria based on parameters defined by functions of the kinetics variables associated with individual tracers would also be of interest. MSE improvements associated with combined analysis of ML-reconstructed data are smaller than those achieved with FBP-reconstructed data. This phenomenon may well have to do with the implicit regularization that ML achieves in low-count frames. Our studies used a common bandwidth for all injections. While this gives some assurance that the resolution of each scan is the same for FBP; it is not clear, given the implicit regularisation achieved by imposition of positivity, if this is necessarily the case for ML. The effect is probably most significant in the repeat studies with many injections. In this setting the total uptake bandwidth selection will lead to very noisy individual scans. The bandwidth selection process ensures voxel-level FBP time-course data are inherently noisier than in ML so the analysis of kinetics for FBP data can end up being worse than the same analysis of ML. The issue is apparent in Table 5. An obvious way to address this issue could be to select bandwidths for individual injections separately. This may allow the noise characteristics of FBP to be more comparable to ML (O'Sullivan;1995;O'Sullivan and Roy Choudhury;2001). The issue merits further investigation - 2021) reports on some initial work in this direction.
For computational reasons we made extensive use of a 1-D scanning model in this work. This allowed comprehensive evaluation of estimators with iteratively reconstructed ML data over a range of count rates and many replicates. While FBP is used very little in practice, we included it in simulations. The relation between 1-D and 2-D FBP results might be used as a vehicle to understand how combined analysis methods of 1-D ML data might perform in a 2-D or 3-D setting. Our simulations here demonstrate that combined analysis will significantly improve the accuracy of NPRM kinetic mapping in multiple injection studies. It will be important to see how these results translate into 3-D and especially for the next generation total-body PET scanners (Badawi et al.;Karp et al.;2020). In this context the emerging techniques for efficient bootstrapping of 3-D multi-frame PET studies (Kucharczak et al.;2020) may facilitate a more direct practical assessment of the benefits of combined analysis in a given patient, without the need for highly sophisticated step-by-step representations of the scanner emission measurements and reconstruction processes.

Appendix:: Derivation of Formula in (11)
Gaussian approximation of estimates of α-coefficients means that as any voxel we can assume α j N α, Σ j where Σ j is a covariance matrix associated with the weighted least squares estimation process. Within the Gaussian framework α j is a sufficient statistic, i.e. it captures the full information about α based on the j'th can data. As a result the optimal (maximumlikelihood) combined estimator of α is a matrix-weighted average of the individual scan estimators, i.e.
where ∑ C = ∑ 1 −1 + … + ∑ J Since λ j ≥ 0, it follows that for any parameter (θ) defined by a linear combination θ = ξ′α i.e. the combined estimator will always have a smaller MSE than the single scan estimate.
More generally suppose that Gaussian approximation applies directly to the kinetic parameter, θ = g(α). This might follow by first-order Taylor series expansion of g(α) about α ( Van der Vaart;2000). With this we would have θ j N θ, σ j 2 wℎere σ j 2 = v′∑ j v (24) and v is the gradient of g evaluated at α. As in (19), θ j is now (to within the approximation) sufficient for estimation of θ based on the information acquired in the j'th scan. The optimal combined scan estimator is a weighted average of these estimates θ C = σ 1 −2 θ 1 + σ 2 −2 θ 2 + … + σ J −2 θ J σ 1 −2 + σ 2 −2 + … + σ J −2 Thus θ C N θ, σ C 2 with σ C 2 = 1 σ 1 −2 + σ 2 −2 + … + σ J −2 . The result in (11) follows from this. Non-parametric residue mapping (NPRM) process.         Relative Performance (6) of NPRM and 2C kinetic mapping in a PET-FDG setting. Based on a 1-D simulation model with source sub-TACs in case I matched to real Breast Cancer data Figure 8 and to their best-fitting 1C model fit in case II. Data reconstruction by both direct (FBP) and iterative (ML) methods.   Voxel-level errors in NPRM kinetics as a function of the local reconstruction (σ i ) and modelling (σ i λ ) errorsc.f. (7), (8). Results based on 1D-PET-FDG simulation.