Propagation of Blood Function Errors to the Estimates of Kinetic Parameters with Dynamic PET

Dynamic PET, in contrast to static PET, can identify temporal variations in the radiotracer concentration. Mathematical modeling of the tissue of interest in dynamic PET can be simplified using compartment models as a linear system where the time activity curve of a specific tissue is the convolution of the tracer concentration in the plasma and the impulse response of the tissue containing kinetic parameters. Since the arterial sampling of blood to acquire the value of tracer concentration is invasive, blind methods to estimate both blood input function and kinetic parameters have recently drawn attention. Several methods have been developed, but the effect of accuracy of the estimated blood function on the estimation of the kinetic parameters is not studied. In this paper, we present a method to compute the error in the kinetic parameter estimates caused by the error in the blood input function. Computer simulations show that analytical expressions we derive are sufficiently close to results obtained from numerical methods. Our findings are important to observe the effect of the blood function on kinetic parameter estimation, but also useful to evaluate various blind methods and observe the dependence of kinetic parameter estimates to certain parts of the blood function.


Introduction
Positron emission tomography (PET) is a functional imaging modality to observe the physiological processes in the body. To conduct a PET scan, positron-emitting radioisotopes, as a tracer, are injected into the living subject (usually into blood circulation). When a positron encounters and annihilates an electron, it emits two gamma rays in reverse directions which will be sensed at two detectors at roughly the same time. Hence it is possible to locate the source along the line of response using a scanner around the subject. The data from the detector is then used to reconstruct an image of the subject [1].
Temporal variation of the tracer concentration can be obtained through dynamic imaging so that the physiological function of the subject can be tracked more accurately. Such a kinetic approach is commonly used in other imaging modalities (e.g., dynamic contrast enhanced MRI [2]) and photodynamic therapy [3]. Dynamic PET, in contrast to static PET, can provide kinetic parameters that are related to physiologic information and is a useful tool for various clinical and research applications [4][5][6][7][8]. A threecompartment model is commonly used in fluoro-2-deoxy-D-glucose (FDG) studies that we use for tumor analysis to simplify the kinetic model of the tracer molecule of interest. In this model, the input C p (t) is the tracer concentration in the plasma, and the output is the time activity curve (TAC). The TACs are obtained by averaging the activity of a known region of interest. Let f (t) denote the TAC of a specific tissue, then the relation between f (t) and the impulse response of a region h(t) can be obtained using diffusion equations.
The differential equations that describe the FDG threecompartment model (Figure 1) are expressed as follows: International Journal of Biomedical Imaging where C p (t) denotes the tracer concentration in the plasma assumed to be spatially constant, C u (t) denotes the concentration of unbounded tracer, C b (t) is the concentration of the bounded tracer, and t denotes the time coordinate. The solution is [8] for t > 0, where " * " denotes convolution and We use the following transformations as commonly done [9] since the transformed parameters are more suitable for our mathematical model: with inverse transforms , Then, for the three-compartment tissue modeling, the relation between f (t) and the impulse response h(t) containing the kinetic parameters is where h(t) is Estimation of the kinetic parameters k 1 , k 2 , k 3 , and k 4 for three-compartment tissue modeling based on f (t) requires that the blood function C p (t) is known. The classical method of arterial sampling to obtain the blood function has several disadvantages: it requires well-trained medical personnel and Figure 1: Three compartments used to model the transfer of the tracer between physical compartments and chemical states.
poses a vital risk to the subject. Therefore, blind methods are developed to estimate the kinetic parameters of the tissue response without knowing input function. In such methods, the input function is estimated along with the kinetic parameters of the tissue impulse model of interest, and thus scale parameters such as k 1 , k 2 , k 3 , and k 4 and any expression using these parameters can only be estimated in a relative sense instead of an absolute one. Several studies have been done in the field, such as maximum likelihood, cross-relation methods, and several others, see [10][11][12][13] and references therein.
Since we have discrete measurements at times t = Δt · n for noisy TAC measurements we can use the following model which can be written as where H (i) is the convolution matrix of the impulse response of the tissue for region i, C p denotes the vector of the blood function, and ε (i) is the noise vector. Stacking different regions of interest together, we can write the equations in the following form We can estimate the kinetic parameters and the blood function by minimizing the following cost function: Several methods for blind kinetic parameter estimation has been proposed, but no study has shown the effect of the errors in the estimated blood on the estimation of kinetic parameters except for our preliminary work [14]. In this paper, we develop a method that can compute this effect in an efficient manner. Our results can be used to calculate the error in the kinetic parameter estimates stemming from the errors in the blood function that is used. Our derivations is based on the implicit function theorem previously used for static PET [15] and the Runge-Kutte methods [16].
Although, these errors in the parameters can be found performing a separate optimization for each of the error combinations in the blood that we want to study, this is a very time-consuming method, considering that the optimization procedure is iterative. This is especially important when we are interested in pixel by pixel kinetic parameter estimation, and/or when the space of the erroneous blood functions we International Journal of Biomedical Imaging 3 want to analyze is large. Based on the results of this paper, the optimization needs to be performed only once, and the error propagation can be calculated very fast based on this single optimization.
Our major contribution in this paper is to construct a mathematical model to derive the error in the kinetic parameter estimates caused by the error in estimation of the blood input function. Our results are conceptually important to observe the effect of the error in the blood function on kinetic parameter estimation, and also practically useful to evaluate various blind methods and observe the dependence of kinetic parameter estimates to certain parts of the blood function such as the peak and tail part. In Section 2, we illustrate the derivation of the mathematical model. Section 3 includes computer simulations that validate the analytical results. Conclusions are presented in Section 4.

Derivation of the Error Propagation for Three-Compartment Tissue Modeling
In this section, we explain how we can calculate the errors in the kinetic parameters for three-compartment tissue modeling due to the error in the blood function. We assume that a unique solution for the blood estimates, at least locally, exists. The estimates of the kinetic parameters a, b, c, d are Our goal is to obtain the errors in the kinetic parameters Δa, Δb, Δc, and Δd stemming from the error in the blood function ΔC p . The first step in this calculation is to calculate the derivatives of the implicit estimator with respect to the elements of C p (t) at a fixed C p (t) point. The second step is to accumulate the error sequentially to derive the ultimate error based on the first step. The first step can be performed by using the chain rule and implicit function theorem [15].
For a solution a, b, c, d that minimizes the cost function the partial derivatives of cost function with respect to the kinetic parameters is zero and three similar equations. Let us define an implicit function which will be convenient for the application of chain rule. We need to use the derivatives of the following equation to derive the expression of the derivatives of the implicit estimator with respect to the elements of C p (t) at a fixed C p (t) point a, b, c, d = g C p = g 1 C p , g 2 C p , g 3 C p , g 4 C p (14) that maps the C p into an estimate [ a, b, c, d]. By applying the chain rule to this equation, we can differentiate the above two equations with respect to C p and obtain and similar equations for b, c, and d can be obtained easily.
Thus we derive the derivative expression that we are interested in, Because of the presence of ∂, b ,c, and d in the terms ∂a/∂ C pn , ∂b/∂ C pn , ∂c/∂ C pn , and ∂d/∂ C pn , we cannot simply integrate (∂a/∂ C pn )Δ C pn , (∂b/∂ C pn )Δ C pn , (∂c/∂ C pn )Δ C pn , and (∂d/∂ C pn )Δ C pn to obtain the errors in the kinetic parameter estimates. However, we can calculate ∂a/∂ C pn , ∂b/∂ C pn , ∂c/∂ C pn , and ∂d/∂ C pn at a fixed point of C pn , and a sequential procedure can be applied to calculate a new value of a, b, c, and d, and we can use them to evaluate ∂a/∂ C pn , ∂b/∂ C pn , ∂c/∂ C pn , and ∂d/∂ C pn at the next C p value until 4 International Journal of Biomedical Imaging the complete range of ΔC p is covered. This procedure can be performed by methods such as Runge-Kutte methods, predictor-corrector method and Richardson extrapolation. , where k i = y i (C pn ), i = 1, 2, 3, 4 for a, b, c, and d, and h is the small step size we define according to the demand on accuracy and speed. This completes the calculation of the errors in a, b, c, and d stemming from the error in estimation of C pn . To summarize the expressions derived in this section provide the errors in the kinetic parameters Δa, Δb, Δc, Δd stemming from errors in the blood ΔC p . The number of steps to derive the ultimate error depends on the step size you specify in different cases to meet the requirement of the accuracy. Then we can use (3) to transform a, b, c, and d to obtain k 1 , k 2 , k 3 , and k 4 . These equations show the computational advantage of the proposed method of calculating the error compared to optimization approach. In the optimization approach, we need to calculate the cost function and the gradients for each of the iterations, whereas here we need to calculate (17) for only a few steps.
The detailed computation of our mathematical model can be seen in the appendix.

Computer Simulations
We apply the proposed mathematical model to . These values are selected to produce a realistic case in [9]. The blood function is produced using the model [12]: with typical values from [12] [  The input blood function can be seen in Figure 2.
We simulate a noisy TAC value as follows: where N is Gaussian Noise with variance 0.01. The noisy TAC model for liver, background, and tumor can be seen in Figure 3. This model adds noise with the power proportional to the activity level to mimic a realistic case. We have performed two sets of experiments. For these three experiments, we set an error margin for C pn and compare the propagated errors in k 1 , k 2 , k 3 , and k 4 of liver and tumor using the derived expressions and optimization. Optimization is a numerical method using the conjugate gradient method (CGM). Table 1: The percentage error in estimated k 1 , k 2 , k 3 , and k 4 of liver affected by a range of erroneous blood functions in three cases.
Error rate of C p Percentage error in parameters 20% 10% First, we test the mathematical model in one dimension and assume that one of the 19 samples of C pn has an error −20% to 20% defined as (C pn (error) − C pn (true))/C pn (true). Figure 4 shows a comparison of the results from the derived expressions and optimization for two arbitrary samples of C p . We observe that the derived expressions provide very accurate approximations of the the errors in k 1 , k 2 , k 3 , and k 4 of liver.
For several of the blind methods, the error in the blood function is not usually confined in a single sample. To illustrate this fact, we have performed a simple simulation where we have estimated the blood function with three different noise realizations. Figure 5 shows that, the initial peak, transition part, and tail section of the estimated blood function is affected differently motivating this type of simulation. Therefore, in our second experiment, we use a blood function with multiple erroneous samples. The blood function is divided into two parts: (i) the initial peak and (ii) the tail part. Based on this grouping three cases of errors are considered: Case I: all samples; Case II: initial peak with samples 1 to 13; Case III: tail part with samples 14 to 19 are erroneous. All erroneous samples have the same error rate ranging from −20% to 20% with the three cases defined before. Figures 6 and 7 show that the results from the derived expressions are very close to ones obtained from numerical optimization.
The following items summarize the effect of the blood function error to the final kinetic parameter estimates.
(a) For Case I, when all 19 samples have the same error rate from −20% to 20%, we can see that the changes in k 2 , k 3 , and k 4 are negligibly small while change in k 1 is relatively large. k 1 of liver drops from 1.0093 to 0.6728, k 1 of background drops from 0.0126 to 0.0084, k 1 of tumor drops from 0.7599 to 0.5067. And k 2 , k 3 , and k 4 data in the Figures 5, 6, and 7 are flat lines. This can be explained with (7) and (8) where scaling in the blood function would not affect c and d but inversely scale a and b.
(b) For Case II, we observe that k 1 and k 3 deviate from the true value considerably, k 4 also deviates from the true value but not as much, while k 2 almost remains the same, indicating that the error in the initial peak affects the estimation of kinetic parameter k 1 , k 3 , and k 4 more than k 2 . In addition to the figures, Table 1 lists the relations between the error rate of C p and the percentage of changes in estimation of the parameters k 1 , k 2 , k 3 , and k 4 of the liver in three cases defined before.
Our conclusion of the relation between the blood function error and the error on the kinetic parameters can be summarized as follows: (a) We can see that the error in the initial peak of the blood input function affects the estimation of kinetic parameter k 1 , k 3 , and k 4 more than k 2 . The parameter k 3 changes 12.5% when error rate of C p becomes 10% while k 1 changes 10% and k 4 changes 4%. (b) And the error in the tail part affects the estimation of kinetic parameter k 3 most, k 4 second but does not affect k 1 and k 2 . The parameter k 3 will change 10% when error rate of C p becomes 10% while k 1 will change 5%.
(c) If the overall blood input function has almost the same error rate, we find that the error in the parameters k 2 , k 3 , and k 4 is negligibly small while k 1 changes relatively large. In this case, only the parameter k 1 will change 10% when error rate of C p becomes 10%. Table 2 lists the relations between the error rate of C p and the percentage of changes in estimation of the parameters k 1 , k 2 , k 3 , and k 4 of the tumor in three cases defined before. Similar conclusion can be drawn from the data of the table.
These simulation results show that the derived expressions provide a very accurate approximation of the errors in the kinetic parameters, and several useful observations related to the effect of the blood function on the kinetic parameter estimates can be made.

Conclusion
Blind identification is recently studied to estimate the kinetic parameters for the compartment without a known blood input function since the arterial sampling of blood input function involves vital risks, requires trained personnel, is not comfortable for the patient in clinical applications, and is difficult to perform in small animals. There are several solutions for blind identification: maximum likelihood methods, cross-relation method, mixture analysis method, and factor analysis of dynamic structures.
Despite several blind methods, how an erroneous blood function affects the final parameter values has not been studied. In this paper, we have derived mathematical expressions that quantify how errors in the blood estimate propagate into errors in the kinetic parameter estimates. Our model is constructed on the base of the implicit function theorem and the Runge-Kutte methods. We first derive the derivative of the kinetic parameters with respect to blood input function at a fixed point. Then we implement the Runge-Kutte approximation to calculate the accumulated errors affected by the gradually increased error in blood input function. The accuracy of the mathematical model can be modified by adjusting the number of steps in the Runge-Kutte approximation as desired.
Computer simulations show that the proposed mathematical model can yield accurate estimates of the errors in k 1 , k 2 , k 3 , and k 4 stemming from the errors in the blood function. We can infer from the results that the error in the initial peak of the blood input function affects the estimation of kinetic parameter k 1 , k 3 , and k 4 more than k 2 . The error in the tail part affects the estimation of kinetic parameter k 3 most, k 4 second but does not affect k 1 and k 2 . If the overall blood input function has almost the same error rate, we find that the error in k 2 , k 3 , and k 4 is negligibly small while k 1 changes relatively large as apparent from equations.
The developed method can quantify the errors in the kinetic parameters for different error combinations in the blood function, without having to perform optimization for each of the error cases to be analyzed. Instead of iteratively optimizing the result of the error using the old method, we can use this analytical method to derive the ultimate error step by step. This would be computationally prohibitive especially for pixel by pixel kinetic parameter estimation, and large ranges of blood error to be analyzed. Future work includes generalization to estimation based on the sinogram instead of reconstructed TAC's and application to real PET data.

Appendix
For a solution a, b, c, d that minimizes the cost function the partial derivatives of cost function with respect to the kinetic parameters is zero Let us define an implicit function which will be convenient for the application of chain rule a, b, c, d = g C p = g 1 C p , g 2 C p , g 3 C p , International Journal of Biomedical Imaging 9 that maps the C p into an estimate [ a, b, c, d]. We can rewrite (13) as By applying the chain rule to this equation, we can differentiate the above two equations with respect to C p and obtain and similar equations for b, c, and d can be obtained easily. We can write these four equations in matrix form Then a simple matrix inversion provides The elements of the matrix and the vector on the right hand side can be calculated based on the cost function and φ is a scalar given by the determinant where h is as defined before. In multidimension, we perform the calculations N times sequentially for every C pn for a single step (A.13)