Patient-adaptive Population-based Modeling of Arterial Input Functions

Kinetic modeling of dynamic PET data requires knowledge of tracer concentration in blood plasma, described by the arterial input function (AIF). Arterial blood sampling is the gold standard for AIF measurement, but is invasive and labour intensive. A number of methods have been proposed to accurately estimate the AIF directly from blood sampling and/or imaging data. Here we consider fitting a patient-adaptive mixture of historical population time course profiles to estimate individual AIFs. Travel time of a tracer atom from the injection site to the right ventricle of the heart is modeled as a realization from a Gamma distribution, and the time this atom spends in circulation before being sampled is represented by a subject-specific linear mixture of population profiles. These functions are estimated from independent population data. Individual AIFs are obtained by projection onto this basis of population profile components. The model incorporates knowledge of injection duration into the fit, allowing for varying injection protocols. Analyses of arterial sampling data from 18F-FDG, 15O-H2O and 18F-FLT clinical studies show that the proposed model can outperform reference techniques. The statistically significant gain achieved by using population data to train the basis components, instead of fitting these from the single individual sampling data, is measured on the FDG cohort. Kinetic analyses of simulated data demonstrate the reliability and potential benefit of this approach in estimating physiological parameters. These results are further supported by numerical simulations that demonstrate convergence and stability of the proposed technique under varying training population sizes and noise levels.


B. Modeling of arterial input data
Image-derived input functions (IDIF) [6]- [13] are obtained by fitting an AIF model to arterial blood pool tracer activity extracted from PET image data, rather than obtained from arterial sampling. Such data are usually extracted from a region of interest (ROI) within a dynamic PET image that captures either a large blood vessel or the left cardiac ventricle, and therefore tend to be noisy. Usually, the average or the highest signals in the ROI are extracted to obtain a time course to characterize the concentration of tracer in arterial blood. IDIF approaches may also require at least one physical blood sample to scale the image-derived AIF curve, and may include tracer uptake information both before and after tracer metabolism. Alternatively, PET signals from multiple regions (not necessarily arterial blood pools) may be leveraged to recover AIF and thereby kinetic parameters, via simultaneous estimation of the input function [14].
information on the patient's population or by averaging fits from the overall cohort) based on individual subject information (depending on availability, historical cohorts have 10-30 blood samples through 60-90 minutes of data, and have a combination of patient weight, body surface area, lean body mass or cerebellar FDG activity) [21], [22]. The choice of time points at which this information is extracted is a likely contributor to variability in estimation of the signal scale. With this approach, the individual is assumed to have the same tracer injection protocol and physiological characteristics as the population, which is not always realistic.
Various mathematical constructions such as Feng's model [23], tri-exponential models [24] and convoluted models [25] have been developed to provide a continuous and noise-free description of the AIF obtained from either clinically sampled (as in typical PBIF strategies) or image-derived (as in typical IDIF strategies) signals, as shown with step 3A of the flowchart in Fig. 1. The tri-exponential construction and its convolved alternative rely on a linear mixture of exponential density functions fit to sample arterial data. The tri-exponential model (TE) fits the AIF shape in two parts: a first linearly increasing piece from the time the AIF starts rising at to its peak, and a second piece modeling the AIF from its peak to the end of study using a linear combination of the three exponential curve. In comparison to Feng's model, the TE model provides more flexibility to describe the AIF. However, neither model accounts for the length of the tracer injection and could lead to poor fit to the initial part of the AIF curve. An adaptation of the TE model was proposed [25] that includes the information of injection duration by convolving the injection profile and linear combination of three exponential curves. This convolved tri-exponential (CTE) model addresses some of the issues of the original TE model, but relies on a boxcar function of fixed duration to describe the injection profile, which used on its own may not represent the pattern of tracer arrival into plasma realistically. This particular limitation has been addressed in some works [26]- [29], by using a convolution of a double Butterworth function with the TE function [28], or of a square wave with a Gamma function [29]. Final individual AIF estimates are obtained by scaling the estimated AIF pattern to the individual physiological level using either image-derived blood level estimates from a reference region (after partial volume correction) or blood sampling data. (Some alternative scaling methods exist that use routine patient information instead; this is discussed further below.) An alternative, compartment model based technique was developed by our group to estimate the AIF from a whole-body tracer circulation model (BCM) [29], [30]. In this representation, the fate of each tracer atom is represented by an eight-state Markov chain process discrete by heart beats, where each state represents one important location in the blood circulation system. To use the BCM, theoretically, time course of the tracer concentration from one of the eight states is required. Since the left ventricle is one of those states represented in the stochastic process, the AIF signal can be approximated by the tracer concentration in the left ventricle observed from PET imaging over time. Otherwise clinically sampled arterial data can be used for evaluation of one of these Markov states. Although high dimensionality of the parameters involved in the model constitutes a computational challenge, this physiologically meaningful representation can yield an accurate AIF estimation.

C. Contribution
Here we consider a modified PBIF-type approach by adjusting a population AIF profile to individual characteristics extracted from either arterial sampling or image-derived data, to achieve reliable patient-adaptive AIF estimation throughout the whole time course. In a sense the proposed population-based projection model (PBPM) combines population profiling (as in a PBIF approach) with individual arterial input data modeling (as in an IDIF approach). As such, it differs from both these strategies, in a number of ways. Unlike IDIF strategies, the PBPM leverages population information, and does not depend on the availability of imaging data (although it can also be applied using such data). Unlike PBIF strategies, it does not average population components into a single individual template pattern, but instead produces a different pattern for each individual, by fitting a linear mixture of population profile components to subject-specific arterial sampling data. Fig.  1 describes the overall process and illustrates how it differs from typical PBIF and IDIF processes. To produce individual fits, the PBPM mixture components are fitted to population data, and a combination of these fits is adapted to each subject, unlike with TE-type methods. This approach could also be adapted for use in step 3B in the flowchart, but no such report was found in the literature. The proposed PBPM allows the use of functions other than exponentials, should a particular tracer require it, and can operate with a different number of components depending on the PET tracer used. This number can be automatically selected; two or three such population components are typically sufficient for individual AIF estimation. The proposed PBPM also provides more flexibility with respect to temporal characteristics of the data. It allows for injection duration to be included in the modeling, which makes it applicable to PET studies with different injection durations. Unlike TE-type constructions, it also does not require specification of two key time points (signal start and peak time), and can adapt to different peak times thanks to its parametrization. A model fitting procedure for the PBPM is also proposed hereafter. A penalized positive constrained least squares method is developed and implemented in R [31] for estimation of individual mixing weights. Statistical regularization is built in to adapt relative weight placed on the individual and population information according to the noise level in a given dataset, to facilitate patient-adaptive fitting of cohort characteristics. This paper focuses on describing aspects related to the PBPM (step 3B in the flowchart), and benchmarking this model against reference methods typically used in step 3A of PBIF and IDIF frameworks. Section II develops the proposed model, and provides details on the methodological developments for fitting the model to clinically sampled arterial data, including regularization and cross-validation for selection of the number of components. Section III illustrates model selection and evaluation, convergence of the procedure and other aspects related to its construction, using both simulated and clinical data comprising of three separate arterial sampling datasets of respectively FDG-, H 2 O-, and FLT-PET imaging studies. Results from these numerical analyses include validation of an open source implementation of the software in R [31] that we made available [32]. Section IV presents a comparative analysis against another three methods of reference, namely the TE, CTE and BCM models. Results on clinical data demonstrate the gain derived by training the proposed AIF representation on population data, as well as its impact on kinetic analysis is evaluated via resampling of uptake templates in one-and two-tissue compartment models.

II. Methods
This section presents the modeling strategy used for patient-specific AIF representation based on PET tracer concentration profiling, including aspects of model fitting and calibration.

A. AIF model specifications
Let the travel time T of a specific atom at a specific arterial sampling site be modeled as the sum of the time (T ir ) for the atom to initially progress from the injection site to the right ventricle (RV) of the heart, and the time (T c ) it spends in circulation, before being sampled at the site: T ir is modeled as a realization from a Gamma distribution (G) with shape and rate parameters α and β, respectively. The circulation time is modeled as a shifted excess circulation time (ξ) with the value of the shift (Δ) being specific to the location of sampled blood-site within the circulatory system If the probability density (f) of ξ is known, it is possible to combine this distribution with the injected dose profile in the right ventricle to evaluate the tracer concentration at any particular time. Under the conventional assumptions of time invariance and linearity of tracer kinetics, a bolus injection at time t = 0 will produce a temporal concentration of tracer atoms in the RV proportional to the distribution G. A more complex temporal injection profile, e.g. specified by an indicator function I(·), would yield an RV concentration profile C RV characterized by an accumulation over the whole injection time, such that C RV (t) = G ⊗ I(t) = ∫ 0 t G(t − s)I(s)ds The evolution of this profile to the concentration profile at the blood site of interest, C P , is governed by the circulation time distribution. This is given by a shifted convolution with the excess circulation travel time density We approximate the distribution f of excess circulation travel time ξ by a linear mixture of J components to obtain an estimate f(ξ) = π 1 e 1 (ξ) + π 2 e 2 (ξ) + ⋯ + π J e J (ξ) where the π 1 , π 2 ,…, π J are mixing fractions of the component densities e 1 , e 2 ,…, e J , with constraint π j > 0, ∀j = 1,…, J. With this representation, an estimate of the concentration profile becomes C P J (t) = ∑ j = 1 J π j ∫ 0 t e j (t − s)G ⊗ I(s − Δ)ds (1) In this expression, the convolution (⊗) of the square wave I() with the Gamma distribution G describes a population injection profile, adjusted to individual subjects by a patient-specific delay Δ. In this approach the population component densities e 1 , e 2 ,…, e J are assumed to be standard exponentials with characteristic rates ϕ 1 > ϕ 2 > ⋯ > ϕ J > 0. While many other possibilities could be considered (e.g. Gamma, Weibull, etc.), a mixture of exponential distributions makes for a simple and flexible model. We presume the mixing components are representative of the entire circulatory system, and functional components e 1 , e 2 ,…, e J and G are held fixed across subjects in the population of study.
Ambiguity may arise in the estimation of the linear coefficients π 1 and π 2 , in cases where the corresponding distribution ranges overlap. To remove this ambiguity, a constraint π 1 > π 2 can be imposed when estimating these linear parameters. This constraint ensures that at early time points, the first component contributes mostly to the peak of AIF representation. A similar constraint may be applied to further mixing coefficients, although in our experience on the data of Section IV, it was not necessary as there was no practical ambiguity between π 2 and π 3 in a three-component mixture applied to dynamic FDG-PET data.

B. Model fitting
Suppose blood time course data from a training population of n subjects {(t ik , z ik , w ik ), i = 1, 2, …, m k , k = 1, 2, …, n} are available, where z ik is the concentration of the ith measurement of the kth subject, measured at time t ik with measurement reliability w ik and m k is the number of time points used in clinical sampling for the kth subject. Each time course z ik is modeled for i = 1, 2, … , m k , k = 1, 2, … , n as where C P (t | Δ k , π k , θ) is defined in the previous subsection with the subject-dependent linear parameter vector π k = c(π 1k , π 2k , … , π Jk ), and subject-independent non-linear parameter vector θ = c(α, β, ϕ 1 , … , ϕ J ). The ε ik 's are assumed to be i.i.d realizations of Gaussian random variables modeling additive noise. We propose to estimate θ by minimizing the objective function defined by l(θ) = ∑ k = 1 n min Δ k , π k ∑ i = 1 m k w ik (z ik − C P (t ik | Δ k , π k , θ)) 2 An algorithm for this optimization is implemented in R [31]. Minimization of l(θ) requires two-step estimation of π k and Δ k , which is described in the following two subsections, respectively.

1)
Estimation of π k given Δ k and θ: The linear mixture is fitted by minimizing, for k = 1 , … , n, criterion with respect to π k given Δ k and θ, where the jth component of the linear mixture for the kth patient, hereafter denoted by A ijk for convenience, can be calculated as where the mixing parameters {π jk } j = 1 J are constrained positive. An R implementation of the Goldfarb-Idnami method [33] is employed to solve this positivity-constrained least-squares problem. We denote the resulting minimizer of l 1 (π k |Δ k , θ) as π k = π(Δ k , θ).
2) Estimation of delay Δ k given θ: Inserting π k (Δ k , θ) in the objective function l 1 , we have the reduced objective function l 2 (Δ k | θ) = l 1 (π k (Δ k , θ) | θ) Given population characteristics θ, a grid search may be carried out to minimize l 2 (Δ k |θ) with respect to Δ k . The measured data may not be coincident with the timing of the injection, and a parameter Δ k , positive or negative, is introduced to account for this. In this study, 21 evenly distributed grid points in the range of [−60 s, +60 s] were used. The resulting minimizer of l 2 is denoted as Δ k .

3) Estimation of population parameter θ:
Finally, the population-specific parameters θ = (α, β, ϕ 1 , … , ϕ J ) are estimated by minimizing Initial values for θ 0 = (α 0 , β 0 , ϕ 1 0 , …, ϕ J 0 ) may be selected based on experience from previous experiments and some physiologic reasoning about the average travel time of each type of tracer in the blood circulation system. Lower and upper bounds θ L = (α L , β L , ϕ 1 L , …, ϕ J L ) and θ H = (α H , β H , ϕ 1 H , …, ϕ J H ), respectively, can be set based on experience from the population of patients receiving the same tracer injection, and physiologic reasoning. Since the scales of the parameters differ considerably, and also to avoid the complexity of constrained minimization, the following parameter transformation may be performed.
Thus, the constrained minimization with respect to θ can be performed as unconstrained minimization of l(ς) with respect to ς. A Gauss-Newton type method was implemented in R to optimize l(ς).

4) Fitting individual curves:
The population parameter estimates θ = arg min θ l(θ) defined in Step (II-B.3) can be used instead of being re-evaluated when fitting time activity curves for "new" patients (i.e. patients outside of the training population of n studies used to obtain θ ), thus reducing the number of parameters required to obtain an individual fit. In such cases, simultaneous estimation of π k and Δ k can be achieved by grid search over G values {Δ k 1 , …, Δ k G }, corresponding estimates {π k 1 , …, π k G } being obtained instantaneously via linear regression.

5) Model regularization:
In cases of low signal-to-noise ratio 1 , regularization of the patient-specific mixing weights π can be applied to control variability of these estimates.
Once the number and parameters of the population components have been determined, a penalized nonlinear least squares approach using the historical linear parameters as a Bayesian penalty for individual fit [34] is used for this approach. Regularization on π = π k parameters for the kth patient is achieved using the following individual objective function: in Section II-B.1. All components in the vectorized mixing weights π are constrained 1 e.g. if using small and/or noisy PET imaging data for AIF estimation to be nonnegative, and π 0 and Ω respectively denote the mean and covariance matrix of linear parameters obtained from population-based estimation. regularization parameter λ controls the trade-off between individual and population data contributions, larger values of λ constraining individual estimates to be closer to the population profile π 0 , which suits higher-noise scenarios. The optimization is a quadratic problem for which we used the Goldfarb-Idnami method [33].
Parameter λ is set by minimizing the generalized cross-validation (GCV) score which is defined as the mean squared error adjusted by the effective degree of freedom: where z ik (λ) is C P (t ik | Δ k , π k , θ) based on minimization of the objective function defined in Eqn. (2), for a given λ.

C. Model selection
Similarly to compartmental analysis of dynamic PET data, the number of components needed in the linear mixture defined in Eqn.
(1) may depend on the tracer used in the studies. More components may be included if the tracer molecules can participate in metabolic activity. For example, in water studies, tracers do not participate in any metabolic activity, and a one-tissue compartmental model is normally used for water dynamic studies. Models using more than one compartment are applied to studies involving other tracers for more complicated circulatory or metabolic systems [4].
Here a leave-one-out approach is developed for selection of the number of components to be used in the proposed linear mixture defined in Eqn. (1) so as to minimize model complexity defined as where J = 2, 3 … is the number of components, C P J ( ⋅ ) defines the predicted AIF based on PBPM with J components given linear and non-linear parameters, and A cross-validated error (Err) may not necessarily be the minimum model complexity Err for the chosen J, but should be reasonably small for all values Err(J) in the range of J.

D. Scale factor estimation
An overarching objective of AIF estimation is to design methods that can rely on imagederived input data instead of invasive sampling. In this context, scaling of PBPM estimates of the AIF obtained from image-derived arterial signals is required to obtain final parametric imaging estimates. This evaluation typically requires at least one blood sample point (from either arterial or venous blood) [15], [35], [36], or an adequately calibrated image-derived measurement from reference tissue (e.g. cerebral tissue) [22], [37], which is not always feasible. Routine individual subject information may alternatively be used to this end. In [37] the authors use FDG uptake measured in cerebellum ROI from a single scan acquisition to normalize a PBIF estimate. In [38] the authors show that body surface area can be used as an alternative to blood sampling in PBIF normalization in brain 11 C-TMSX PET studies.
In [39], the authors use injection dose and weight information to scale PBIF estimates in 11 C-DPA-713 PET studies. Here we propose a regression model for the scale factor that does not require arterial blood sampling or a large blood-pool in the field of view. This model uses injection dose X ID , blood volume X BV (estimated from non-invasive patient information on body weight, height and total blood volume [40]) and tail height of the AIF curve (π 3 , the linear coefficient for the third component in the PBPM) as independent variables: log(S) = β 0 + β 1 X ID + β 2 X BV + β 3 π 3 + noise (4)

III. Evaluation
This section introduces the three clinical datasets used for evaluation of the proposed PBPM model. Model selection, evaluation of the impact of training population size on model fitting, and scaling of the final AIF estimates are covered in these analyses. Complementary simulations based on the sampling data also illustrate convergence of the procedure at varying noise levels.

A. Data
The proposed model was applied to three separate datasets, comprising blood sampling data obtained during the acquisition of respectively 105 FDG-PET studies, 39 H 2 O-PET studies, and 32 FLT-PET studies, all acquired at the University of Washington [9], [29], [41]- [44]. The data were from both healthy subjects and patients with cancer conditions. All FLT patients had cancer (brain, breast, lung, sarcoma). Of the 105 FDG blood curves used in the analysis, 12 were from normal human subjects with the rest from patients with cancer, including brain, lung, liver and colon cancer patients. Eighteen of the 39 water curves came from normal subjects. Durations of FDG injections were either 1 or 2 minutes (58 and 33 cases respectively), bolus injection for H 2 O was assumed to be 5 seconds, and duration of FLT injection was 1 minute. The datasets are summarised in Table I All water studies were manual bolus injections. Only the initial 7 FDG studies were manual injections over 1 min. The other FDG studies had either a 1 min or 2 min programmed infusion using a syringe pump. All FLT injections were a 1 min injection from an infusion pump. A saline flush was used after every injection. Ten studies from each of these 3 sets were used for PBPM fitting of the population components θ. The rest of the data were used to obtain individual fits from each of the models used in the following comparative analyses. The normalized AIFs are presented in Fig. 2. Each curve typically includes three parts: a plateau of 0 concentration at the begin, followed by a rapid increase after injection, and an exponential decay. No imaging data were used in the following analyses.

B. Model selection
The clinical datasets were used to estimate the population-based parameters in the proposed PBPM model for varying numbers of components. Figure 3 depicts the cross-validated error Err defined from Err(J), for the FDG-and H 2 O-PET cohorts respectively, after carrying out the model selection procedure described in Section II-C. For the FDG set an elbow was clearly identified at 3 components, indicating a 3-component mixture to be optimal in the sense of loss Err(J). Three exponentials with different rates convolved with a Gamma distribution were thus estimated from the population data to approximate the travel time distribution of the FDG and FLT tracers, yielding a reasonable representation in PBPM (which aligns with the traditional two-tissue compartment modeling of FDG molecules in the blood circulation system). For the H 2 O set, there was no clear elbow and a 2-component mixture seemed the optimal choice according to this selection method. Since, unlike FDG, H 2 O is not involved in the metabolic process, it seems reasonable that the number of components used to describe H 2 O AIF might be one less than FDG. Note that this model selection rule is only indicative, and using an over-specified mixture (i.e. using more components) may also mechanically impact the bias-variance tradeoff achieved by the model.

C. Model fitting
Model fitting was analysed on both clinical data (Section III-C.1) and simulated data created using the clinical data as templates (Section III-C.2), to illustrate the procedure and demonstrate its convergence and reliability at varying noise levels.

2) Model validation and assessment:
A set of complementary simulation-based analyses were conducted to validate the model and its implementation. A first analysis was carried out on AIF templates generated by adding random noise to the PBPM model curve, to validate the implementation and evaluate model performance at varying noise levels. Details on the simulation settings used for these analyses are provided in Appendix A. Figure 5 presents one example simulated curve and its fit for each noise level. Figure 6 presents the performance of PBPM on 400 simulated AIFs where medium noise level is as above when the noise pattern is set to be similar to the residuals obtained from fitting the PBPM model to the directly sampled AIFs. The RMSE distributions showed overall improvement in model fit as noise level decreases, as expected. The parameter estimation error distributions in Fig. 6 showed an increase in estimation efficiency as noise level decreases. The figure also illustrates that model performance using nonlinear parameters estimates for θ = (α, β, ϕ 1 , ϕ 2 , ϕ 3 ) was comparable to that obtained when using the true population parameter values. These results demonstrate that the PBPM approach is consistent and accurate, and validate its implementation [32].

D. Sensitivity to training population size
As for the training population size required to fit the mixture components, behaviour of the PBPM model may vary with the PET tracer used. An experiment on the FDG cohort was carried out, using random subsets ranging from n = 5 to n = 60 studies taken from the overall set as training population studies, and monitoring cross-validated model fit errors. The results from this analysis, shown in Fig. 7, indicated that model fitting was relatively insensitive to this parameter, with comparable error distributions obtained when using population parameter estimated across the range of population sizes.

E. Scale factor estimation
Results on 91 FDG AIF studies from the scale factor estimation procedure described in Section II-D are depicted in Fig. 8 (4 patients were removed from the set of 95 curves used in other analyses, due to missing information). They show significant alignment between estimated scale factors (Eqn. (4)) and actual scale factor values obtained directly from arterial sampling data (Pearson correlation ρ = 0.83, two-sided test with p < 1e −10 under the null hypothesis H 0 : ρ = 0).

IV. Results
This section provides output of comparative analyses between the proposed PBPM and reference approaches, using the three clinical datasets introduced in Section III-A. The benefit of using basis components fitted to separate population data is also evaluated hereafter.

A. Comparative analysis
PBPM-derived AIF estimates were compared to alternatives obtained from another three models described in Section I-B, namely a whole-body blood circulation model (BCM) [29], [30], a straightforward tri-exponential model (TE) [5], [23], and a convolved tri-exponential model (CTE) [25]. For each tracer, the PBPM was fitted using a random sample of 10 arterial blood sampling curves for estimation of the population parameters θ = (α, β, ϕ 1 , …, ϕ J ) (and not for any other step), and the rest of the data for estimation of individual AIF signals, i.e. of {(π k , Δ k ), k = 1, … , n} (and not for any other step). In other words, the 10 randomly selected "population curves" were used as "historical data" in this framework. These were only used for the purpose of fitting the PBPM, and not the BCM, TE and CTE models. A comparative analysis between these four models was carried out in terms of the sum of mean squared errors obtained for each of the n scaled AIF fits {z i, k = C P (t i, k | π k )} k = 1 n , defined by This fit error was cross-validated in the comparative analysis, using (5) for comparison, where z i, k * = η * (t i, k | π k * ) and η * (t i, k | π k * ) is the set of optimal parameter estimates for the data {(t j, k , z j, k )} j = 1, j ≠ i m .
The boxplots of cross-validated error distributions of Fig. 9

B. Gain from a population-based approach
An initial numerical experiment was carried out on simulated data to assess the potential gain offered by employing a population-based strategy over a fully individual approach to estimate the basis components when fitting the PBPM to sampling data. The same simulation process defined by Eqn. (11) was applied to generate 210 simulated curves from the PBPM templates. A random set of 10 of these curves were used as a training population to obtain estimates θ of the nonlinear parameters θ. The remaining 200 curves were fitted once following the proposed methodology (using this estimate θ to generate PBPM fits), and independently a second time by re-estimating the nonlinear parameters θ individually (i.e. without using the estimate derived from the training population). We refer to the former AIF estimate as the PBPM fit, and the latter AIF estimate as the individual fit. Fit errors were calculated for each of the 200 curves and cross-validated using leave-one-out cross-validation. Figure 10 shows the distributions of cross-validated fit errors obtained from each of the PBPM and individual fits, at varying noise levels. The results show a clear gain resulting from the population-based methodology.
A second analysis was carried out on the clinical data to confirm these findings. Figure 11 illustrates the gain obtained by using a population-based strategy to estimate the nonlinear parameters θ for the PBPM basis component profiles, over a fully individual fitting process where all model parameters (π, Δ, θ) would be estimated from the single patient sampling data. A one-sided paired Wilcoxon test, under the null hypothesis that the average crossvalidated PBPM error is not less than the average cross-validated individual fit error, confirmed the gain achieved when using estimates θ obtained from a separate training population (p < 0.0001).

C. Impact on kinetic modeling
AIF estimation is of interest in the context where there is calculation of kinetic parameters. Impact of the proposed PBPM approach on parametric imaging was evaluated in that sense via Monte Carlo simulation of H 2 O-and FDG-PET time activity curves generated from one-tissue compartment (1C) and two-tissue compartment (2C) models respectively [4], [45]- [48], where  [50] for modeling of grey matter activity in cerebral MRglu studies). This process yielded a set of respectively 39 and 91 unique template uptake curves. These templates were normalised to 1 so estimation of a scale factor was not required in this experiment. Finally, simulated uptake curves were obtained for each scenario by adding Gaussian noise with standard deviation σ = ϕ z t , true with ϕ = 0.001 for H 2 O and ϕ = 0.04 for FDG, to 10 replicates of each of these template curves. This yielded respectively 390 and 910 simulated time activity curves for these two settings.
Estimates of flow (K 1 ), rate k 2 , and distribution volume V D = K 1 /k 2 for the H 2 O scenario, and of flux (K i = K 1 k 3 /(k 2 + k 3 )), distribution volume (V D = K i /k 3 ) and rates K 1 , k 2 and k 3 for the FDG scenario, were compared across the range of 1C and 2C model fits obtained using AIF estimates from BCM, TE, CTE and PBPM for the H 2 O-and FDG-PET cohorts. Note that in the FDG case here, quantity V D [51], [52] differs from the distribution volume V t = K 1 /k 2 (1 + k 3 /k 4 ) of total ligand uptake in tissue relative to total concentration of ligand in plasma typically of interest for 2-TC models, since we assume k 4 = 0 [53]. For FDG, the distribution of non-displaceable compartment relative to total concentration of ligand in plasma V nd = K 1 /k 2 is also reported on. (In the case of water, these quantities are equivalent.) Delay was held constant for the purpose of this analysis, and the compartment models were optimized with respect to (K 1 , k 2 , k 3 , V B ) (with k 3 = k 4 = 0 in the 1C model, and k 4 = 0 in the 2C model) via nonlinear estimation. Figure 12 shows the distributions of percentage errors in estimating these parameters using varying AIF estimates, and Table III provides the p-values of corresponding one-tailed Wilcoxon (Mann-Whitney) tests to assess bias from these distributions, under the null hypothesis of a strictly lower bias for the alternative model, compared to PBPM (i.e. H 0 : Err PBPM ≥ Err other , after adjustment of error signs where appropriate for the comparison). In complement, "mirror" tests assessing H 0 : Err PBPM ≤ Err other indicated that PBPM bias was significantly greater than BCM bias for estimation of V nd and rate k 3 in the FDG setting (p < 0.005 for both). No other such tests yielded significance. In particular, TE and CTE did not yield significantly lower bias compared to PBPM in any scenario.
In summary, the PBPM model is a viable option for evaluation of physiological parameters in the FDG and H 2 O settings, where it led to either statistically comparable or significantly reduced bias in every scenario for all kinetic parameters, with the only exception of BCM for estimation of V nd and k 3 in FDG. These results illustrate its competitiveness and reliability for kinetic analysis across different tracers, against the most effective methods of reference.

V. Discussion
This paper introduced a novel population-based projection model (PBPM) that combines strategies used in both image-derived input function (IDIF) and population-based input function (PBIF) approaches in order to achieve flexible and reliable AIF estimation throughout the whole time course, by fitting a mixture of population-informed uptake patterns to individual patient imaging data. In this sense this methodology can be seen as the projection of individual patient PET tracer kinetics characteristics onto a basis of population uptake Using a few population-based components, the proposed method can incorporate knowledge of the injection duration into the model fit, thus not requiring the injection protocol to be the same in all PET studies, unlike common PBIF methods.
A comparative analysis was carried out that included two reference models, namely the tri-exponential (TE) and convolved tri-exponential (TCE) models, and a whole-body blood circulation model (BCM) as references for FDG-, H 2 O-and FLT-PET datasets, as well as simulated datasets obtained from resampling of the clinical information. Results demonstrated that the proposed PBPM has overall competitive or lower AIF estimation error and greater adaptability to different tracers compared to the other three methods, and illustrated its consistency at varying noise levels. Application of the proposed approach to parametric imaging, using one-and two-tissue compartment kinetic modeling scenarios, was also evaluated. In this experiment, the PBPM method overall yielded appropriate characterisation of kinetic rates as well as flow, flux and distribution volume, and yielded either comparable or significantly improved bias over the other models (except for two cases where alternate BCM yielded statistically significant improvement over PBPM), showing improved adaptivity to different tracers.
A natural end-goal is to apply AIF modeling to image-derived data, either because clinical sampling is not available or to facilitate non-invasive evaluation. As part of the technical contributions of this paper, regularization of the model should confer robustness against increased noise levels typically found in imaged blood-pool data. A region of interest from PET imaging of the left ventricle could be used in this sense, as the contrast of concentration of tracer from ventricles and their surrounding muscle is usually distinct. Follow-on work is planned around future access to datasets that combine sampling and imaging data to carry out such an evaluation. This will also yield the opportunity to evaluate the proposed PBPM method compared to IDIF techniques used commonly with major PET tracers.
Where only individual venous blood samples and no reliable arterial blood information are available, the proposed technique may still be considered, for example after conversion of the venous signals into arterial scales. Some venous input function (VIF) approaches were proposed that apply an arteriovenous transform scheme to estimate AIF, for example in 11 C-labelled studies [54].
The requirement for scale factor estimation, discussed in Section II-D, impacts all modeling techniques, including the proposed PBPM strategy, but it is not yet clear how. Analysing this impact will be the scope of follow-on work.
A critical minimum number of PET studies required to use the proposed PBPM technique with future studies remains to be determined and could vary with the PET radiotracer used. Numerical experiments demonstrated robustness of this model in this aspect, with stable model fit errors for even very few population studies.
This study has not discriminated the data in terms of healthy and cancer subjects, nor in terms of disease types. Although we do not anticipate significant differences in AIF profiles across such strata, using imaged data may introduce differences in PBPM estimates of AIF profiles. Further analyses would be required to explore this.

VI. Conclusion
We proposed a statistical model of the arterial input function of a patient that leverages arterial uptake profile information at both the cohort and subject levels. Results from simulations demonstrated numerical stability and consistency of the approach at varying noise levels. Improvement in accuracy of the AIF representation over reference AIF models was assessed on three clinical cohorts, each imaged with one of three routinely used PET tracers (namely, FDG, H 2 O and FLT). The proposed model either outperformed or remained competitive against three reference techniques, in terms of AIF estimation RMSE, for these tracers. Experiments further demonstrated that this representation of the AIF is viable for kinetic analysis of physiological parameters such as flow, flux and distribution volume from one-and two-tissue compartment models. Overall the proposed model yielded statistically comparable or significantly lower bias against the three alternate models for determination of kinetics from typical FDG template data in almost all comparisons. The proposed methodology can be incorporated in either population-based or image-derived input function estimation frameworks; the latter will be the scope of follow-on work assessing this model for non-invasive AIF estimation.

A. Validation study on model-derived template
In the first analysis of Section III-C.2, three mixture components were fit to the FDG-PET imaging cohort to obtain estimates of the Gamma and exponential parameters with parameters set as follows: α = 2.1784, β = 7.8741, ϕ 1 = 88.3992, ϕ 2 = 0.1488, ϕ 3 = 0.0088 .
Step 3 of each method is an AIF model fitting step. Blood sampling for step 6A is not used in some IDIF approaches (hence the dashed line); when used, fewer points are typically required compared to blood sampling used in step 6B. Using individual image-derived data (1A) at step 3B(ii) of the PBPM approach (instead of arterial sampling data) would require an extra step in order to scale the output to a physiological level, as in step 5 of an IDIF method. Normalised kinetic PET time activity curves acquired at the University of Washington [9], [29], [41], [43], [44]:     Top: comparison of cross-validated errors obtained when using a training population strategy to estimate the nonlinear parameters θ that determine the structure of the PBPM basis components (left), over a fully individual approach where the parametric vector θ is estimated directly from the same individual sampling data (right), showing a lower error tends to be achieved from a population-based strategy (one-sided paired Wilcoxon test p < 0.0001). Bottom: direct comparison based on percent differences between these PBPM and individual fit errors ((E indiv − E PBPM ) ∕ E PBPM × 100).