Exact parameter identification in PET pharmacokinetic modeling using the irreversible two tissue compartment model

Abstract Objective. In quantitative dynamic positron emission tomography (PET), time series of images, reflecting the tissue response to the arterial tracer supply, are reconstructed. This response is described by kinetic parameters, which are commonly determined on basis of the tracer concentration in tissue and the arterial input function. In clinical routine the latter is estimated by arterial blood sampling and analysis, which is a challenging process and thus, attempted to be derived directly from reconstructed PET images. However, a mathematical analysis about the necessity of measurements of the common arterial whole blood activity concentration, and the concentration of free non-metabolized tracer in the arterial plasma, for a successful kinetic parameter identification does not exist. Here we aim to address this problem mathematically. Approach. We consider the identification problem in simultaneous pharmacokinetic modeling of multiple regions of interests of dynamic PET data using the irreversible two-tissue compartment model analytically. In addition to this consideration, the situation of noisy measurements is addressed using Tikhonov regularization. Furthermore, numerical simulations with a regularization approach are carried out to illustrate the analytical results in a synthetic application example. Main results. We provide mathematical proofs showing that, under reasonable assumptions, all metabolic tissue parameters can be uniquely identified without requiring additional blood samples to measure the arterial input function. A connection to noisy measurement data is made via a consistency result, showing that exact reconstruction of the ground-truth tissue parameters is stably maintained in the vanishing noise limit. Furthermore, our numerical experiments suggest that an approximate reconstruction of kinetic parameters according to our analytic results is also possible in practice for moderate noise levels. Significance. The analytical result, which holds in the idealized, noiseless scenario, suggests that for irreversible tracers, fully quantitative dynamic PET imaging is in principle possible without costly arterial blood sampling and metabolite analysis.


Introduction
Positron emission tomography (PET) is a non-invasive clinical technique that images the four dimensional spatio temporal distribution of a radio tracer in-vivo.In quantitative dynamic PET imaging, after tracer injection, several 3D PET images at different time points are acquired and reconstructed.The injected tracer is supplied to the tissues via the arteries and capillaries.As a response to this "arterial tracer input", the tracer is exchanged with tissues.This exchange can include reversible and/or irreversible binding and eventually metabolization of the tracer.
Using the right tracer, the time series of reconstructed PET images reflecting the tissue response, a measurement of the arterial tracer input, and a dedicated pharmacokinetic model, it is possible to generate images reflecting certain physiological parameters.Depending on the tracer, these parametric images can reflect e.g.blood flow, blood volume, glucose metabolism or neuro receptor dynamics.
Pharmacokinetic modeling in PET is commonly performed using compartment models, where the compartments usually reflect tissue subspaces.Example for such subspaces are the extra cellular spaces where the tracer is free or bound. 1 The dynamics between the arterial blood and tissue compartments is typically described using ordinary-differential-equation (ODE) models.For PET tracers with irreversible binding, such as [ 18 F]Fluorodeoxyglucose (FDG) [14] or [ 11 C]Clorgyline [10], the irreversible two tissue compartment model can be used to describe the tracer dynamics, see Figure 1 for a scheme of this model.
The identification of the kinetic parameters describing the tissue response to the arterial tracer supply in a given region is commonly done using the following input data: 1.The tracer concentration in tissue C tis (t).This quantity, which is equal to the sum of all the tracer concentrations in all extra-vascular compartments, can be directly obtained from the time series of reconstruction PET images -either on a region-of-interest (ROI) or at voxel level.
2. The arterial "input function" or plasma concentration C art (t) of the non-metabolized tracer describing the concentration of tracer in the arterial blood plasma that is supplied to tissue and available for exchange (and metabolism).A direct measurement of C art (t) is complicated.Typically, it is based on an external measurement of a time series of arterial blood samples taken from a patient, e.g. using a well counter.Unfortunately, the total arterial blood tracer concentration C bl (t) obtained from the well counter measurements usually overestimates C art (t) since the measured activity of the blood samples also includes activity from radioactive molecules that are not available for exchange with tissue because of (i) parts of the radio tracer being bound to plasma proteins and (ii) activity originating from metabolized tracer molecules that were transfered back from tissue into blood.Measuring the contributions of the latter two effects, summarized in the parent plasma fraction f (t) = C art (t)/C bl (t), requires further advanced chemical processing and analysis of the blood samples and is thus time consuming and expensive.
To avoid the necessity of arterial blood sampling, which itself is a very challenging process in clinical routine, many attempts have been made to derive the arterial input function directly from the reconstructed PET images (also called "image-based arterial input function") by analyzing the tracer concentration in regions of interest of the PET images containing arterial blood, such as the left ventricle, the aorta, or the carotic arteries.Note, however, that by using any imagebased approach for the estimation of the arterial input function, the contributions of tracer bound to plasma proteins and metabolized tracers cannot be determined.In other words, any image-based approach can only estimate C bl (t) instead of C art (t) since no measurement of the parent plasma fraction f (t) is available.Motivated by the problem, we consider the question whether the identification of kinetic parameters using data from multiple anatomical regions of interest and the irreversible two tissue compartment model is possible without having access to the common parent plasma fraction f (t) and/or the common total arterial blood tracer concentration C bl (t).
In existing literature on modeling approaches for quantitative PET, see for instance [17] and [16] for a review, this question has been addressed from a computational perspective.The work [17], for instance, accounts for a low number of measurements of the arterial concentration of non-metabolized PET tracer by using a non-linear mixed effect model for the parent plasma fraction, i.e., the parameters defining the parent plasma fraction are modeled as being partially patient-specific and partially the same for a population sample.Moreover, the general idea of jointly modeling the tissue response in different anatomical regions to obtain unknown common 2 The irreversible two tissue compartment model The irreversible two tissue compartment model describes the interdependence of the concentration of a radio tracer in the arterial blood plasma and in the extra-vascular compartment, where the latter is further decomposed in a free and a bound compartment.Note that in the irreversible model, once the radio tracer has reached the bound compartment, it is trapped.A visualization of the model is provided in Figure 1.
We denote by C art : [0, ∞) → [0, ∞) the arterial plasma concentration of the non-metabolized PET tracer.Further, for any anatomical region i = 1, . . ., n, we denote by the free and the bound compartment of the tracer in region i, respectively, and by C i tis = C i fr + C i bd we denote the sum of the two compartments in region i.Using the irreversible two tissue compartment model, the interaction of these quantities is described by the following system of ordinary differential equations (ODEs): Here, the parameters K i 1 , k i 2 and k i 3 are the tracer kinetic parameters that define the interaction of the different compartments in region i.
Our goal is to identify the parameters K i 1 , k i 2 and k i 3 for each i = 1, . . ., n.For this, we can use measurements of the C i tis (t l ) at different time-points t 1 , . . ., t T obtained from reconstructed PET images at different time points after tracer injection.Further, the parameter identification typically relies on additional measurements related to C art .Here, the standard procedure is to take arterial blood samples and to measure the total activity concentration of the arterial blood samples, given as C bl : [0, ∞) → [0, ∞).The relation of the total concentration C bl to the arterial plasma concentration C art of non-metabolized tracer is described via an unknown parent plasma fraction function As described above, to obtain f (t) and thus C art (t), a time-consuming and costly plasma separation and metabolite analysis of the blood samples has to be performed.
In order to realize the parameter identification for the ODE model (S) in practice, the involved functionals need to be discretized, e.g., via a suitable parametrization.For the arterial concentration C art , it is standard to use a parametrization via polyexponential functions as defined in the following.
Definition 1 (Polyexponential functions).We call a function g polyexponential of degree p ∈ N if there exist λ i , µ i ∈ R for 1 ≤ i ≤ p where the (µ i ) p i=1 are pairwise distinct and λ i ̸ = 0 for all i such that We write deg (g) = p and call the zero-function polyexponential of degree zero.By P p we denote the set of polyexponential functions of degree less or equal to p, and by P = ∪ ∞ p=0 P p the set of polyexponential functions (of any degree).
Remark 2. It obviously holds that P is a vector space, and even a subalgebra of C ∞ (R).It is also worth noting that, as direct consequence of the Stone-Weierstrass Theorem, polyexponential functions are dense in the set of continuous functions on compact domains.Thus, they are a reasonable approximation class also from the analytic perspective.
Modeling C art as polyexponential function already defines a parametrization of the resulting solutions of the ODE system (S).As we will see in Lemma 5 below, the following notion of generalized polyexponential functions is the appropriate notion to describe such solutions.
Definition 3 (Generalized Polyexponential class).We call a function g generalized polyexponential if it is of the form where the P 1 , . . ., P l ̸ ≡ 0 are polynomials of degree m 1 − 1, . . ., m l − 1, respectively, and the µ j are pairwise distinct constants.We denote the class of such generalized polyexponential functions with polynomials of degree at most m 1 − 1, . . ., m l − 1 by We define the degree of g by In case m 1 = • • • = m l = 1 we write P deg(g) for the resulting polyexponential class.
The next two results, which follow from standard ODE theory, provide explicitly the solutions (C fr , C bd ) of the ODE system (S), once in the general case and once in case C art is modeled as polyexponential function.
Lemma 4. Let C art : [0, ∞) → [0, ∞) be continuous, and let the parameters K i 1 , k i 2 and k i 3 be fixed for i = 1, . . ., n such that k i 2 + k i 3 ̸ = 0.Then, for each i = 1, . . ., n, the ODE system (S) admits a unique solution (C i fr , C i bd ) that is defined on all of [0, ∞), and such that C i tis = C i fr + C i bd is given as Proof.Fix i ∈ {1, . . ., n}.From the equation for C i fr in (S) it immediately follows that This in turn implies that and, consequently, that Lemma 5.If C art ∈ P p is given as λ j e µj t , then, for i ∈ {1, . . ., n}, C i tis of Lemma 4 is given as Proof.This follows immediately by inserting the representation of C art in (1).
Remark 6 (Sign of exponents µ j ).Note that for the ground truth arterial concentration C art , we will always have µ j < 0 (in particular µ j ̸ = 0) for all j = 1, . . ., p, since otherwise this would imply the unphysiological situation that C art (t) → c ̸ = 0 for t → ∞.

Unique identifiability
In view of Lemma 5 from the previous section, it is clear that the question of unique identifiability of the parameters K i 1 , k i 2 and k i 3 from measurements of C i tis (t l ) at time points t 1 , . . ., t T is related to the question of uniqueness of interpolations with generalized polyexponential functions.A first, existing result in that direction is as follows.
Lemma 7 (Roots of generalized polyexponential functions).Let P 1 , . . ., P l be polynomials of degree m 1 − 1, . . ., m l − 1 such that at least one of them is not identically zero, and let the constants µ 1 , . . ., µ l be pairwise distinct.Then the function As a consequence of the previous proposition, we now obtain the following unique interpolation result for generalized polyexponential functions.
Lemma 8 (Unique interpolation).Let m 1 , . . ., m p , T ∈ N be such that Then, for any choice of tuples Based on this result, we now address the question of uniquely identifying the parameters of the ODE system (S) from time-discrete measurements C i tis (t 1 ), . . ., C i tis (t T ) with i = 1, . . ., n and measurements C bl (s 1 ), . . ., C bl (s q ).For this, we first introduce the following notation.
Definition 9 (Parameter configuration).We call the parameters p, n ∈ N, ((λ j , µ j )) λ j e µj t a configuration of the irreversible two tissue compartment model if λ j ̸ = 0 for j = 1, . . ., p, the µ j , j = 1, . . ., p are pairwise distinct, and, for i = 1, . . ., n, C i tis = C i fr + C i bd with (C i fr , C i bd ) the solution of the ODE system (S) with arterial concentration C art and parameters K i 1 , k i 2 , k i 3 .Central for our unique identifiability result will be the following technical assumption on a parameter configuration (p, n, ((λ j , µ j )) p j=1 , (( For any j 0 , there are at least three regions i s , s = 1, . . ., 3, where As the following lemma shows, this assumption holds in case our measurement setup comprises sufficiently many regions where the parameters k i 3 and k i 2 + k i 3 are pairwise distinct.This is reasonable to assume in practice, and also it is a condition which is to be expected: Our unique identifiability result will require a sufficient amount of different regions, and different regions with the same tissue parameter do not provide any additional information on the dynamics of the ODE model.Lemma 10.Assume that there are at least p + 3 regions i 1 , . . ., i p+3 , with p ≥ 1, where each the k is 3 and the k is 2 + k is 3 are pairwise distinct for s = 1, . . ., p + 3. Then Assumption (A) holds.
Proof.For z ∈ R, note that is a polynomial in z of degree at most p − 1. Hence it can admit at most p − 1 distinct roots.Now since there are at least p + 3 regions where each the k is 3 and the k is 2 + k is 3 are pairwise distinct, for at least four of them, say i 1 , . . ., i 4 , z = k is 2 + k is 3 cannot be a root of the above polynomial.Further, for those four regions, since the k is 3 are pairwise distinct, for any given µ j0 , at most one can be such that µ j0 + k is 3 = 0.As a consequence, the remaining three are such that the conditions of Assumption (A) hold.
Based on Assumption (A), we now obtain the following proposition, which is the technical basis for our subsequent results on unique identifiability.
Proof.Take (p, n, ((λ j , µ j )) p j=1 , (( , Cart ) to be two configurations as stated in the proposition, such that in particular for l = 1, . . ., T .Now since µ j ̸ = 0 for all j = 1, . . ., p, we obtain the following representation of C i tis : In particular, for any region i ∈ {1, . . ., n}, the coefficients of e µj t for j = 1, . . ., p in this representation are given as either The latter can only be zero if µ j + k i0 3 = 0, which can happen for at most one ĵ by the (µ j ) j being pairwise distinct.Since p ≥ 3 by assumption, this implies in particular that C i tis is a non-zero function for any i.As a consequence of ( 5), the condition T ≥ 2(p + 3) and the unique interpolation result of Lemma 8, this implies that Ci tis is a non-zero function, such that in particular Ki 1 ̸ = 0 for all i ∈ {1, . . ., n}.Together with the assumption that ki 3 ̸ = 0 for all i, we also obtain that μj ̸ = 0 for all j, since otherwise Ci tis would have a non-zero coefficient of t.As a consequence, also Ci tis admits a representation as in (6).Uniqueness of the exponents (µ j ) p j=1 .As first step, we now aim to show that p = p (in particular λj ̸ = 0 for all j) and that (up to re-indexing) µ j = μj for all j = 1, . . ., p.
We start with a region i 0 ∈ {1, . . ., n}.In this region, as argued above, the coefficients of the e µj t for j = 1, . . ., p can be zero for at most one ĵ.Since further at most one j 0 can be such that µ j0 = −( ki0 2 + ki0 3 ), the unique interpolation result of Lemma 8 applied to C i0 tis and Ci0 tis yields that p ≥ p − 2 ≥ 1 and that (up to re-indexing) µ j = μj for all j / ∈ { ĵ, j 0 }.Now as a consequence of Assumption A, we can pick a region i This means that the coefficients of both e µ ĵ t and e µj 0 t in the representation of C i1 tis as in (3) are non-zero.Again by the µ j being pairwise distinct, this implies that p ≥ p − 1 and that (up to re-indexing) either µ ĵ = μĵ or µ j0 = μj0 .Case I. Assume that µ ĵ = μĵ .If also µ j0 = μj0 (and p ≥ p) we are done with this step, so assume the contrary.Now as a consequence of Assumption A, we can pick i 2 , i 3 and i 4 to be regions where µ j0 +k i l 3 ̸ = 0 and either µ j0 3 )t in the representation of C i l tis is non-zero.The fact that µ j0 + k i l 3 ̸ = 0, together with µ j0 ̸ = μj0 by assumption, further yields that µ j0 = −( ki l 2 + ki l 3 ) for all l = 2, 3, 4. Now we argue that in each region i l with l = 2, 3, 4, it must hold that either To this aim, we make another case distinction for a fixed l ∈ {2, 3, 4}.
Case I.A Assume that there exists j l with k i l 2 +k i l 3 +µ j l = 0. From the fact that this can happen at most for one j l and that λ j l ̸ = 0, it follows that the coefficient of te −(k i l 2 +k i l 3 )t is non-zero.Consequently, it follows from the unique interpolation result that k i l 2 + k i l 3 = ki l 2 + ki l 3 as claimed.Case I.B Assume that k i l 2 + k i l 3 + µ j ̸ = 0 for all j.This means that μj = µ j ̸ = −(k i l 2 + k i l 3 ) for all j ̸ = j 0 .But since the coefficient of e −(k i l 2 +k i l 3 )t is non-zero, by the unique interpolation result it must hence hold that either −(k for all l = 2, 3, 4, one of the two cases must happen at least twice.By uniqueness of the k i l 2 + k i l 3 for l = 2, 3, 4, only 3 can happen twice.On the other hand, since µ j0 = −( ki l 2 + ki l 3 ) for all l = 2, 3, 4, this yields that at least two k i l 2 + k i l 3 coincide, which is a contradiction.Hence Case I is complete.Case II.Assume that µ j0 = μj0 .In this case, interchanging the role of µ j0 and µ ĵ , we can argue that µ ĵ = μĵ exactly as in Case I.
Uniqueness of at least three of the exponents k i 2 + k i 3 .Let i 0 be any region such that either µ j0 + k i0 2 + k i0 3 = 0 for some j 0 ∈ {1, . . ., p} (such that the coefficient of te −(k i l tis is non-zero, and note that, according to Assumption A, at least three such regions exist.Case I. Assume that there exists j 0 ∈ {1, . . ., p} such that k i0 2 + k i0 3 + µ j0 = 0.This implies that the coefficient of te −(k i 0 2 +k i 0 3 )t is non-zero and, consequently, already that 3 ) must match some exponent in the representation of Ci0 tis .It cannot match any of the μj = µ j since k i0 2 + k i0 3 + µ j ̸ = 0 for all j = 1, . . ., p, hence again k i0 2 + k i0 3 = ki0 2 + ki0 3 follows.Uniqueness of at least three of the exponents k i 2 , k i 3 .First note that for any i ∈ {1, . . ., n} where ki 2 + ki 3 = k i 2 + k i 3 , from the unique interpolation result, it follows that for all j = 1, . . ., p. Indeed, in case k i 2 + k i 3 + µ j = 0, it follows from the coefficients of te −(k i 3 ) as claimed.In the other case, the equality (7) follows directly from the coefficients of e µj t in C i tis and Ci tis being equal.Now let i 0 be any region where ki0 , and for which we want to show that k i0 2 = ki0 2 and k i0 3 = ki0 3 .Again we consider several cases.Case I. Assume that there exists j 0 ∈ {1, . . ., p} such that µ j0 + k i0 3 = 0.In this case, it follows from (7) that also µ j0 + ki0 3 = 0 (note that λj0 ̸ = 0 and Ki0 1 ̸ = 0 since p = p), hence k i0 3 = ki0 3 and, consequently, k i0 2 = ki0 2 holds.Case II.Assume that µ j + k i0 3 ̸ = 0 for all j.In this case, using Assumption A and the previous step, we can select i 1 to be a second region where again ki1 We have two cases.
Case II.k3 + µ j ̸ = 0 for all j = 1, 2, 3.In this case we can reformulate (12) to obtain for all j = 1, 2, 3.In particular, this yields ), this implies µ 1 = µ 2 and hence a contradiction.Thus, using that z 1 (µ 2 + k3 ) ̸ = z 2 (µ 1 + k3 ) we we can reformulate the previous equation to obtain Now, in equation ( 12) for j = 3, replacing ζ K1 by the equality ( 14) for j = 2 and plugging in the expression ( 15) for k2 we obtain, after some reformulations, k3 Using the definition of the z j in (12) we derive that the factor after k3 in ( 16) corresponds to the term which is nonzero by the µ j being pairwise distinct.Thus, again plugging in the definition of the z j in ( 12) and rearranging the terms in ( 16) yields k3 = k 3 after some computations and, consequently, also k2 = k 2 and ζ K1 = K 1 by the previous considerations.
As a consequence, the remaining ζ Ki 1 , ki 2 , ki 3 considered in this final part of the proof are uniquely determined as ζ Ki The previous result shows that, already under knowledge of C i tis (t l ) for i = 1, . . ., n and sufficiently many distinct time-points t l , the coefficients k i 2 , k i 3 and the coefficients K i 1 can be determined uniquely and uniquely up to a constant, respectively.Considering the ODE system (S), it is clear that this result cannot be improved in the sense that the constant factor of K i 1 cannot be determined without any knowledge of C art (since one can always divide all K i 1 by a constant and multiply C art by the same constant).
In case one aims to determine all parameters of a given configuration uniquely, some additional measurements related to C art are necessary.It is easy to see that a single, non-zero measurement of C art , for instance, would suffice.Indeed, given the value of a ground truth Ĉart (ŝ) ̸ = 0 at some time-point ŝ, the equality C art (ŝ) = Ĉart (ŝ) = Cart (ŝ) together with the result from Proposition 11 immediately imply that ζ = 1 such that all parameters are uniquely defined.
In current practice, indeed measurements of C art are obtained via an expensive blood-sample analysis, and used for parameter identification, see for instance [17].As discussed in the introduction, however, in contrast to obtaining measurements of C art , it is much simpler to obtain measurements of the total concentration C bl , where C art = f C bl with the unknown parent plasma fraction f .
As the following result shows, measurements of C bl only are indeed sufficient to uniquely identify all remaining parameters, provided that one has sufficiently many measurements in relation to a parametrization of f .To formulate this, we need a notion of parametrization of the parent plasma fractions.
Definition 12 (Parametrized function class for parent plasma fractions).For any q ∈ N, we say that a set of functions F q ⊂ {f : R → R} is a degree-q parametrized set if for any f, f ∈ F q and λ ∈ R it holds that λf − f attaining zero at q distinct points implies that λ = 1 and f = f .Simple examples of degree-q parametrized sets of functions are polynomials of degree q − 1 that satisfy f (x 0 ) = c for some given x 0 , c ∈ R with c ̸ = 0 or polyexponential functions of degree q/2 (if q is even) that satisfy f (x 0 ) = c for some given x 0 , c ∈ R with c ̸ = 0.The latter is a frequently used type of parametrization for parent plasma functions (where f (0) = 1 is required), see for instance [17].
Proposition 13.In the situation of Proposition 11, assume in addition that f, f : R → R are parent plasma fractions contained in the same degree-q parametrized set of functions, and are such that with s 1 , s 2 , . . ., s q being q different points, and C bl (s l ) ̸ = 0 given for l = 1, . . ., q.Then, all assertions of Proposition 11 hold with ζ = 1, and further Proof.Proposition 11 already implies that Cart = ζC art .Using that, by assumption, we obtain (ζf − f )(s l ) = 0 for l = 1, . . ., q.Since f, f : R → R are parent plasma fractions contained in the same degree-q parametrized set, this implies that ζ = 1 and f = f as claimed.
The following theorem now summarizes results of the previous two propositions in view of practical applications.Theorem 14.Let (p, n, ((λ j , µ j )) p j=1 , ((K i 1 , k i 2 , k i 3 )) n i=1 , (C i tis ) n i=1 , C art ) be a ground-truth configuration of the irreversible two tissue compartment model such that 1. p ≥ 3, n ≥ 3 and K i 1 , k i 2 , k i 3 > 0 for all i = 1, . . ., n, 2. There are at least p + 3 regions i 1 , . . ., i p+3 where each the k is 3 and the k is 2 + k is 3 are pairwise distinct for s = 1, . . ., p + 3.
If further f : [0, ∞) → [0, ∞) is a ground-truth parent plasma fraction in a degree-q parametrized set of functions and f : [0, ∞) → [0, ∞) is a parent plasma fraction in the same degree-q parametrized set of functions such that C art (s l ) = f (s l )C bl (s l ) and Cart (s l ) = f (s l )C bl (s l ) for l = 1, . . ., q, with the s 1 , . . ., s q pairwise distinct and C bl (s l ) ̸ = 0 given, then ζ = 1 and Proof.This is an immediate consequence of Lemma 10 and Proposition 11: Indeed, Lemma 10 ensures that the assumptions of Proposition 11 are satisfied provided that 1.) and 2.) hold.In case p ≤ p the result immediately follows from Propositions 11 and 13.In case p > p it follows from interchanging the roles of the two configurations and again applying Propositions 11 and 13.
Remark 15 (Interpretation for practical application).Besides putting some basis assumptions on the ground truth-configuration and requiring positivity of the metabolic parameters, the previous theorem can be read as follows: If one obtains a configuration that matches the measured data, it can be guaranteed to coincide with ground-truth configuration if at least p + 3 of the found terms ki 2 + ki 3 and ki 3 are pairwise distinct.Remark 16 (Generalization for nontrivial fractional blood volume).Following standard approaches in quantitative PET modeling, we assume here that the PET images provide exactly the tissue concentration C tis .A more realistic model would be that the voxel measurements provide a convex combination of the tissue and blood tracer concentration given by , where fbv with 0 ≤ fbv ≤ 0.05 describes the fractional blood volume.In case the parameter fbv is known and C bl is available at the same time points as the PET image measurements, our results cover also this setting.The general case, where both fbv and C bl are unavailable, can be addressed by similar techniques as in the proof of Proposition 11.Here, the idea would be to employ a polyexponential parametrization also for C bl , and assuming enough measurements of C PET to be available in order to apply the unique interpolation result of Lemma 8.One would further have to ensure positivity of C bl , the intial condition C bl (0) = 0 and conditions on the attenuation f = C art /C bl such as monotonicity and limiting conditions with respect to time approaching zero and infinity, respectively.These requirements imply corresponding conditions on the parameters of C art and C bl .

A Tikhonov approach for parameter identification with noisy data
In the previous section we have established that, under appropriate conditions, the parameters of the irreversible two tissue compartment model in regions i = 1, . . ., n can be obtained uniquely from measurements C tis (t i ), i = 1, . . ., T and measurements C bl (s i ) for i = 1, . . ., q.While this result shows that parameter identification is possible in principle, it considers the idealized scenario of exact measurements.In this section, we consider the situation of noisy measurements, for which we develop a Tikhonov approach for stable parameter identification.As main analytic results of this section, we show i) a stability result, i.e., that the proposed Tikhonov approach is stable with respect to (noise) variations in the measurements (see Theorem 18) and ii) a consistency result, i.e., that in the limit of vanishing noise, solutions of the Tikhonov approach converge to the ground truth parameters (see Theorem 20).
As first step, we define a forward model that maps the unknown parameters to the available measurement data.To this aim, we define the arterial concentration as mapping Further, we define a parametrized parent plasma fraction as mapping where M ⊂ R q is some (finite dimensional) parameter space and F q is a degree-q parametrized set of functions.
Remark 17 (Parent plasma fraction example).A classical model for the parent plasma fraction (see [16] for different models), that we will also use in our numerical experiments below, is the biexponential model Here 2 and the degree of F q is q = 4.
In addition to the parameters of the functions modeling the arterial concentration and the parent plasma fraction, the forward model also includes the parameters 3 ) for i = 1, . . ., n compartments.With this, the unknown parameters are summarized by (λ, µ, m, K i , . . ., K n ) and we denote by X = R p × R p × R q × R 3×n the resulting parameter space with norm where ∥ • ∥ 2 denotes the Euclidean norm.Given measurement points t 1 , . . ., t T for C tis and s 1 , . . ., s q for the total concentration C bl , those parameters are mapped forward to a measurement space Y = R n×T +q , again equipped with the Euclidean norm ∥ • ∥ Y = ∥ • ∥ 2 via the function where, for x = (λ, µ, m, K 1 , . . ., K n ) and Here C tis C art (λ, µ) , K i , denotes the solution of the irreversible two tissue compartment ODE model (S) with parameters K i and arterial concentration C art (λ, µ).Note that F 2 depends on the data C bl that must be obtained from blood samples or PET measurements, which we assume to be given throughout this section.A further adaption of the model to include also C bl as possibly noise measurement is possible with the same techniques as below, but will be omitted for the sake of simplicity.Now denoting by Ĉi tis (t 1 ), . . ., Ĉi tis (t T ) for i = 1, . . ., n measurements corresponding to the ground-truth parameters, our goal is to find parameters (λ, µ, m, K 1 , . . ., K n ) such that Accounting for the fact that the given parameters are perturbed by measurement noise, i.e., we are actually given (C i tis ) δ (t l ) with we address the parameter identification problem via a minimization problem of the form min Here 0 ∈ R q is a q-dimensional vector of zeros, C δ tis summarizes the available measurements for C δ tis , i.e., and λ, μ, m, K1 , . . ., Kn is an initial guess on the ground truth parameters.The above approach corresponds to Nonlinear Tikhonov-Regularisation, for which stability and consistency results can be ensured as follows.
Theorem 18 (Well-posedness and stability).Let the parent plasma fractions f m be such that the mapping m → f m (t) is continuous for any t ∈ [0, ∞).Then, for any given datum C δ tis , the minimization problem (21) admits a solution.Moreover, solutions are stable in the sense that, if (C δ k tis ) k is a sequence of data converging to some datum C δ tis , then, any sequence of solutions (x k ) k of (21) with data (C δ k tis ) k admits a convergent subsequence, and the limit of any convergent subsequence x is a solution of (21) with data C δ tis .
Proof.Since X and Y are finite dimensional and D(F ) is obviously closed, this follows from classical results in regularization theory, see for instance [4, Theorem 10.2] provided that F is continuous.
We start with continuity of F 1 as in (19), the first component of F .For this, it suffices to show that the mapping from the the parameter (λ, µ, K 1 , . . ., K n ) to C i tis (t), with t ∈ [0, ∞) fixed, is continuous, which, in turn, follows from the representation of C i tis (t) as in (1) if, for any g ∈ L 2 (0, t) and any sequence (λ l , µ l ) l converging to (λ, µ) it holds that By Hölder's inequality, the latter follows from C art λ l , µ l → C art (λ, µ) in L 2 (0, T max ), which, in turn, follows via the Lebesgue dominated convergence theorem from point-wise convergence of C art λ l , µ l and the fact that |C art λ l , µ l | on [0, t] can easily be bounded by a constant independent of l.
Regarding F 2 as in (20), the second component of F , continuity immediately follows from continuity of (λ, µ) → C art (λ, µ)(t) and m → f m (t) for any t ∈ [0, ∞) fixed, where the latter holds by assumption.
Remark 19 (Continuity of m → f m (t)).Note that the assumption of continuity of m → f m (t) is only necessary since we allow for arbitrarily parametrized parent plasma fractions; it holds in particular for the biexponential model of Remark 17 and will typically hold for any reasonable parametrization.
At last in this section we now establish a consistency result.
Theorem 20 (Consistency).Let (p, n, (( λj , μj )) p j=1 , (( Ki Ĉart ) be a groundtruth configuration of the irreversible two tissue compartment model satisfying the assumptions of Theorem 14, and let f m ∈ M be a ground-truth parent plasma fraction.
With x = ( λ, μ, m, K1 , . . ., Kn ) the corresponding parameters and ŷ := F (x) the corresponding measurement data, let y δ k be any sequence of noisy data such that ∥ŷ −y Then, any sequence of solutions (x k ) k of (21) with data y δ = y δ k and α = α k such that α k → 0 and δ 2 k /α k → 0 as k → 0 admits a convergent subsequence.Any limit x = (λ, µ, m, K 1 , . . ., K n ) of such a subsequence, such that the corresponding parameter configuration satisfies the assumptions of Theorem 14, coincides with x.Further, if any limit of a convergent subsequence corresponds to a parameter configuration satisfying the assumptions of Theorem 14, then the entire sequence Proof.This is a consequence of Theorem 14, which ensures that there is a unique x ∈ X with F (x) = ŷ, and of classical results from regularization theory, see for instance [4,Theorem 10.3].
Remark 21 (Interpretation of the consistency result).When choosing p ≥ 3 and n ≥ 3, and given the definition of D(F ) as in (18), the above consistency result together with the unique reconstructability result of Theorem 14 can be interpreted as follows: Whenever the parameters ) n i=1 corresponding to a limit x of (x k ) k are such that at least p + 3 of the parameters k i 3 and k i 2 + k i 3 are pairwise distinct, then one can ensure that x = x.
Remark 22 (Multi-parameter regularization).The setting of (21) and the subsequent results on well-posedness and consistency can be generalized to incorporating different regularization parameters for the different norms and data terms, see for instance [6], which is reasonable given the fact that the parameters might live on different scales, and given the fact that the noise level of different measurements over time might be different.
Remark 23 (Model Variations).Currently, in the setting of (21), the parameters 3 ) n i=1 are bounded away from zero by ϵ > 0. For the (µ j ) p j=1 , we currently do not pose any constraints even though, as mentioned in Remark 6, only the choice µ j < 0 is reasonable from a physiological perspective.Likewise, C art as parametrized in (17), does not necessarily satisfy C art (0) = 0.These two conditions, however, can be easily incorporated in the model via the additional constraint µ j ≤ −ε for some ε ≥ 0 and via setting λ p = − p−1 j=1 λ j , respectively.

Numerical solution algorithm
In this section, we provide proof-of-concept numerical experiments that illustrate the analytic unique identifiability results of Sections 3 and 4 also numerically.

Experimental setup
We consider the following experimental setup: As ground truth data, we consider n = 3 different anatomical regions, where, based on the realistic values provided in [8, Table 1], the parameters are chosen as Here, region 1 corresponds to the frontal cortex, region 2 to the temporal cortex and region 3 to the occipital cortex in the human brain.We model the true arterial concentration by the triexponential function given by and the parent plasma fraction for the biexponential model by for 0 ≤ t ≤ T max where T max = 3750 sec = 62.5 min.See Figure 2 for the visualisation of the parent plasma fraction, the arterial concentration and the measured total and additive concentrations.We assume that we conduct a PET examination for T = 25 time frames where, equidistantly distributed, six take place in the first minutes, four in the following two minutes, two in the next two minutes, three in the next 7.5 minutes and finally ten in the remaining 50 minutes.Regarding C bl , we assume that measurements C bl (s 1 ), . . ., C bl (s q ) are available at the same timepoints 0 = t 1 , . . ., t T at which measurements of C tis are available, i.e., q = T .This reflects the situation that an estimate of the total arterial tracer concentration is obtained from the PET images (rather than blood sampling), a technique that for which many recent works exist [18].In view of Propositions 11 and 13, the above experimental setting satisfies the assumptions such that unique identifiability from noiseless data can be guaranteed.We summarize the unknown parameters of this setting by where x † denotes the ground-truth parameters as specified above.For a given number of measurements T ∈ N, the data is summarized in vectorized form via where, abusing notation, F : R 18 → R 3T +T is a vectorized version of the forward model of ( 18).For the subsequent numerical experiments, we will also consider noisy versions of the data, denoted by y δy with δ y being the noise level.Those are defined by adding Gaussian noise with zero mean and variance δ 2 y 3T to each of the first 3T entries of y † (note that no noise is added to the zero-entries), i.e., (y δ ) l = y † l + η l with η l ∼ N 0, for l = 1, . . ., 3T , such that where E(•) denotes the expectation.
As we are dealing with locally convergent methods, it will also be important to choose a reasonable initial guess x 0 for the algorithm.In order to test the performance of the algorithm in dependence on how close the initial guess is to the true solution, we employ the following steps to obtain perturbed initial guesses.Given a level of perturbation δ x , we define where σ ∼ Unif ({−1, 1}), i.e., is uniformly distributed on {−1, 1} and γ ∼ N δ x , 1 4 δ x , i.e., is Gaussian distributed with mean δ x and variance δ x /4.This results in a expected squared deviation of x 0 from x † as by where σ i ∼ Unif ({−1, 1}) and γ i ∼ N δ x , 1 4 δ x for i = 1, . . ., 18 are independent random variables.

Algorithmic implementation
In order to numerically solve the non-linear parameter identification, we employ the iteratively regularized Gauss-Newton method of [1], see also [4,Section 11.2].This is a standard method for solving non-linear inverse problems.It is related to the Tikhonov approach discussed in Section 4 in the sense that similar results on stability and convergence/consistency (under appropriate source conditions) can be obtained, see for instance [2,5], but different to the Tikhonov approach, regularization is achieved by early stopping of the algorithm rather than adding an additional penalty term to the data-fidelity term.Early stopping has the advantage that, using an estimate of the noise level of the data, the discrepancy principle [4, Section 4.3] can be used to determine the appropriate amount of regularization, without requiring multiple solutions of a minimization problem as would be the case with the Tikhonov approach.
Given an initial guess x 0 ∈ D(F ) and a sequence of regularization parameters (α k ) k such that where c α > 1 is some constant, the iteration steps of the iteratively regularized Gauss-Newton method for k = 0, 1, 2, . . .are given as where Those iteration steps are repeated until the discrepancy principle is satisfied, that is, until ∥F x δ k − y δ ∥ Y ≤ τ δ holds for the first time, where δ is an estimate of ∥y δ − y † ∥ Y and τ > 1 with τ ≈ 1 is a parameter.The iterate x k is then returned as the approximate solution of F (x) ≈ y δ .
Remark 24 (Guaranteed convergence).Since the parameter identification problem addressed here is highly non-linear, global convergence guarantees for any numerical solution algorithm are out of reach.For the iteratively regularized Gauss-Newton method together with the discrepancy principle, as considered here, at least local convergence guarantees can be obtained as long as a particular source condition, i.e., a regularity condition on the ground truth solution holds, see [5] for details.
In a practical application, the iteration (27) is combined with a projection on D(F ), which is a closed, convex set for which the projection is explicit (we denote the projection map by P D(F ) ), see [9,Theorem 4] for corresponding results on convergence of such a projected method.Together with this, we arrive at the algorithm for solving F (x) ≈ y δ as provided in Algorithm 1, where we set ϵ = 10 −3 for defining D For the regularization parameters (α i ) i we choose the ansatz for i ∈ N where a = 800 and b = 1/5 are fixed parameters.Besides fulfilling the decay conditions of (26), this choice is motivated by the goal of penalizing deviations from the initial guess rather strongly at early iterations (a large), and avoiding an exploding condition number of the matrix , that needs to be inverted at each iteration, during later iterations (b rather small).
For the realization of the forward operator F and the adjoint of its Fréchet-Differential the main idea is to vectorise the computations and omit expensive for-loops.For that, one may exploit that the entries of F (x) and F ′ [x] * , which mostly consist of integral type entities, may be computed analytically.The elementary components are of the form e µs ds, e (k2+k3)(s−t) e µs ds, se (k2+k3)(s−t) e µs ds and se µs ds.The latter may be computed by hand applying integration by parts.The corresponding terms, which depend on a combination of  Experiments where carried out for each combination of δ y and δ x as above.In order to obtain representative results, each experiment was carried out 100 times.Those experiments where the algorithm diverged (i.e., no improvement compared to the initialization was achieved) were dropped (see Table 1 for the number of dropped experiments for each parameter combination) and, among the remaining ones, the one whose performance was closest do the median performance of all repetitions was selected for the figures below.
Figures 4 to 7 show the result obtained for different choices of δ y = 0, 10 −4 , 10 −3 , 10 −2 .In those figures, the top and bottom rows depict the evaluation of the residual value ∥F (x k ) − y δ ∥ Y and the relative error ∥x k −x † ∥ X ∥x † ∥ X over iterations, respectively, both with a logarithmic scale for the vertical axis.The different columns show results for the different choices of δ x = 0.15, 0.1, 0.05, 0.01.The values of ρ opt and ρ d (the latter only for δ y > 0), together with the respective iterations, are provided at the top of plots of the second row.Note that, while the algorithm was always run until a fixed, maximal number of iterations (300 for δ y = 0 and 200 for δ y = 0) for obtaining the figures, in practice the iterations would be stopped by the discrepancy principle for δ y > 0. For δ y > 0, the red lines always indicate the levels of the residual value (top row) and iteration number (bottom row) where the algorithm would have stopped according to the discrepancy principle.
In Figure 4, which considers the noiseless case, one can observe that the ground truth pa-Figure 4: Performance of IRGNM for δ y = 0 and different initial guesses rameters x † are approximated very well for all levels of δ x .This confirms our analytic unique identifiability result also in practice.
Considering Figures 5 to 7, one can observe that, in the low noise regime, a good approximation of the ground truth x † is still possible across different values of δ x .For higher noise, a good initialization (i.e., a small value of δ x ) is of increasing importance and, in some cases, in particular for δ y = 10 −2 and the parameters of the parent plasma fraction f , the ground truth can not be recovered reasonably well.
In a second set of experiments, we consider the situation that not only C bl , but also measurements of the values of C art are available at the time-points t 1 , . . ., t T .While it is possible in practice to obtain those values via blood sample analysis, this procedure is time consuming and expensive, such that it is a relevant question to what extent such samples improve the identifiability of the tissue parameters.
Again, for each combination of δ y and δ x , each experiment was carried out 100 times, the number of divergent experiments is shown in Table 1 and  While for δ y ∈ {0, 10 −4 }, both versions yield acceptable results, for δ y ∈ {10 −3 , 10 −2 } knowledge of C art enables a good approximation of the ground truth parameters in situations where this was not possible with knowing only C bl .This indicates that, as one would expect, the problem of also identifying f from measurements is significantly more difficult and there is a benefit (at least for this solution methods) in measuring C art via blood samples.

Conclusions
In this work, we have shown that most tissue parameters of the irreversible two-tissue compartment model in quantitative PET imaging can, in principle, be recovered from standard PET measurements only.Furthermore, a full recovery of all parameters is possible provided that suf-    ficiently many measurements of the total arterial concentration are available.This is important, since it shows that parameter recovery is possible via using only quantities that are easily obtainable in practice, either directly from the acquired PET images or with a relatively simple analysis of blood samples.While these results consider the idealized scenario of noiseless measurements, we have further shown that standard Tikhonov regularization applied to this setting yields a stable solution method that is capable of exact parameter identification in the vanishing noise limit.
These findings open the door to a comprehensive numerical investigation of parameter identification based on only PET measurements and estimates of the total arterial tracer concentration, using real measurement data and advanced numerical algorithms.While the numerical results in this paper provide a first indication that it can be possible to transfer our analytical results to concrete applications, we expect that a comprehensive effort is necessary to obtain a practically usable numerical framework for parameter identification.This will be the next step of our research in that direction.

Figure 1 :
Figure 1: Irreversible two tissue compartment model.The boxes around C art , C bd , C fr and C tis represent the respective concentrations and the arrows between them the directional exchange with the respective rate depicted above or below the corresponding arrow.The box around the C bd and C fr indicates the extra-vascular measurements associated with C tis .

Figure 2 :
Figure 2: The evolutions in the first two plots display the ground true parent plasma fraction and arterial concentration.The third and fourth plot correspond to the resulting simulated measurements of the total and additive concentrations depending on the parent plasma fraction and the arterial concentration.

Figure 3 :
Figure 3: Visualization of magnitudes occurring in different time-intervals (i.e., columns) of the data y †

Table 1 :
Number (full setting/setting with known C art ) of experiments (out of ten) not included in the final evaluation due to divergence.δ x = 0.15 δ x = 0.1 δ x = 0.05 δ x = 0 the Figures show the experiment whose performance was closest to the median performance of the non-divergent experiments.Results are shown in Figures 8 to 11, where the quantities shown in and above the plots are the same as in Figures 4 to 7. It can be observed that the performance with known C art is improved compared to the situation where only C bl is known, across all choices of δ y and δ x .

Figure 8 :
Figure 8: Performance of IRGNM for δ y = 0 and different initial guesses and known parent plasma fraction f

Figure 9 :Figure 10 :
Figure 9: Performance of IRGNM for δ y = 10 −4 and different initial guesses and reduced operator

Figure 11 :
Figure 11: Performance of IRGNM for δ y = 10 −2 and different initial guesses and reduced operator