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

Authors: Yafang Cheng; İmam Şamil Yetik
Journal: International Journal of Biomedical Imaging
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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–8]. A three-compartment 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 Cp(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 three-compartment model (Figure1) are expressed as follows:  (1)dCu(t)dt=k1Cp(t)-(k2+k3)Cu(t)+k4Cb(t),dCb(t)dt=k3Cu(t)-k4Cb(t),
where Cp(t) denotes the tracer concentration in the plasma assumed to be spatially constant, Cu(t) denotes the concentration of unbounded tracer, Cb(t) is the concentration of the bounded tracer, and t denotes the time coordinate. The solution is [8]  (2)Cu(t)=k1c-d[(k4-d)e-dt+(c-k4)e-ct]*Cp(t),Cb(t)=k1k3c-d[e-dt-e-ct]*Cp(t),
for t>0, where “*” denotes convolution and  (3)d=(k2+k3+k4)2-(k2+k3+k4)2-4k2k42,c=(k2+k3+k4)2+(k2+k3+k4)2-4k2k42.
We use the following transformations as commonly done [9] since the transformed parameters are more suitable for our mathematical model:
(4)a=k12Δ(k2-k3-k4+Δ),b=k12Δ(-k2+k3+k4+Δ),c=0.5(k2+k3+k4+Δ),d=0.5(k2+k3+k4-Δ),Δ=|(k2+k3+k4)2-4k2k4|
with inverse transforms
(5)k1=a+b,k2=ac+bda+b,k3=ab(c-d)2(a+b)(ac+bd),k4=cd(a+b)ac+bd.
Then, for the three-compartment tissue modeling, the relation between f(t) and the impulse response h(t) containing the kinetic parameters is
(6)f(t)=h(t)*Cp(t),
where h(t) is
(7)h(t)=ae-ct+be-dt.Figure 1
Three compartments used to model the transfer of the tracer between physical compartments and chemical states.Estimation of the kinetic parametersk1, k2,  k3, and k4 for three-compartment tissue modeling based on f(t) requires that the blood function Cp(t) is known. The classical method of arterial sampling to obtain the blood function has several disadvantages: it requires well-trained medical personnel and 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 k1, k2, k3, and k4 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–13] and references therein.Since we have discrete measurements at timest=Δt·n for noisy TAC measurements we can use the following model  (8)f(i)(n)=h(i)(n)*Cp(n)+ε(i)(n),
which can be written as
(9)f⃗(i)=H(i)Cp⃗+ε⃗(i),
where H(i) is the convolution matrix of the impulse response of the tissue for region i, Cp 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  (10)F⃗=HCp⃗+ε⃗.We can estimate the kinetic parameters and the blood function by minimizing the following cost function:(11)R=‖F⃗-ĤĈp⃗‖2.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 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 Section2, we illustrate the derivation of the mathematical model. Section 3 includes computer simulations that validate the analytical results. Conclusions are presented in Section 4.


## 2. 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 parametersâ, b̂, ĉ, d̂ are
(12)[â,b̂,ĉ,d̂]=arg⁡min⁡⁡â,b̂,ĉ,d̂R([â,b̂,ĉ,d̂],Ĉp).
Our goal is to obtain the errors in the kinetic parameters Δa, Δb, Δc, and Δd stemming from the error in the blood function ΔCp. The first step in this calculation is to calculate the derivatives of the implicit estimator with respect to the elements of Cp(t) at a fixed Cp(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 solutiona,b,c,d that minimizes the cost function the partial derivatives of cost function with respect to the kinetic parameters is zero  (13)0=∂∂aR([â,b̂,ĉ,d̂],Ĉp)|a=â
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 Cp(t) at a fixed Cp(t) point
(14)[â,b̂,ĉ,d̂]=g(Ĉp)=[g1(Ĉp),g2(Ĉp),g3(Ĉp),g4(Ĉp)]′
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
(15)0=∂2∂a2R(g(Ĉp),Ĉp)∂∂Ĉpg1(Ĉp)+∂2∂a∂bR(g(Ĉp),Ĉp)∂∂Ĉpg2(Ĉp)+∂2∂a∂cR(g(Ĉp),Ĉp)∂∂Ĉpg3(Ĉp)+∂2∂a∂dR(g(Ĉp),Ĉp)∂∂Ĉpg4(Ĉp)+∂2∂a∂ĈpR(g(Ĉp),Ĉp),
and similar equations for b, c, and d can be obtained easily.Thus we derive the derivative expression that we are interested in,(16)[∂g1(Ĉp)∂Ĉpn∂g2(Ĉp)∂Ĉpn∂g3(Ĉp)∂Ĉpn∂g4(Ĉp)∂Ĉpn]=[∂2∂a2R∂2∂a∂bR∂2∂a∂cR∂2∂a∂dR∂2∂b∂aR∂2∂b2R∂2∂b∂cR∂2∂b∂dR∂2∂c∂aR∂2∂c∂bR∂2∂c2R∂2∂c∂dR∂2∂d∂aR∂2∂d∂bR∂2∂d∂cR∂2∂d2R]-1×[-∂2∂a∂ĈpnR-∂2∂b∂ĈpnR-∂2∂c∂ĈpnR-∂2∂d∂ĈpnR].
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 Cp value until the complete range of ΔCp is covered. This procedure can be performed by methods such as Runge-Kutte methods, predictor-corrector method and Richardson extrapolation.(17)s1,i=hui(Cpn,a,b,c,d),s2,i=hui(Cpn+h2,a+s1,12,b+s1,22,c+s1,32,d+s1,42),s3,i=hui(Cpn+h2,a+s2,12,b+s2,22,c+s2,32,d+s2,42),s4,i=hui(Cpn+h,a+s3,1,b+s3,2,c+s3,3,d+s3,4),yi(Cpn+h)=yi(Cpn)+s1,i6+s2,i3+s3,i3+s4,i6,
where ki=yi(Cpn), 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 ina, b, c, and d stemming from the error in estimation of Cpn. 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 ΔCp. 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 k1, k2, k3, and k4. 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.


## 3. Computer Simulations

We apply the proposed mathematical model to simulated dynamic PET data witht=[0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0,2.5,3.0,3.5,7.0,10.0,15.0,20.0,30.0,60.0,90.0,120.0] min, where the ending time of each frame is used, and three sets of kinetic parameters for liver, background, and tumor: the first set of liver: k1=0.8, k2=0.98, k3=0.012, k4=0.017 (a=0.79, b=0.01, c=0.9922, d=0.0168); the second set of background: k1=0.01, k2=0.05, k3=0.001, k4=0.0001 (a=0.0098, b=0.00019683, c=0.0510, d=0.000098035); the third set of tumor: k1=0.598, k2=0.680, k3=0.05, k4=0.001 (a=0.5569, b=0.0411, c=0.7301, d=0.00093142). These values are selected to produce a realistic case in [9]. The blood function is produced using the model [12]:(18)Cp(t)=(b1t-b2-b3)eλ1t+b2eλ2t+b3eλ3t
with typical values from [12](19)[b1,b2,b3,λ1,λ2,λ3]=[900,22,21,-4.1,-0.12,-0.01]min⁡-1.
The input blood function can be seen in Figure 2.Figure 2
The input blood function.We simulate a noisy TAC value as follows:(20)Fnoisy=Foriginal⋅(1+N),
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.Figure 3
The observed noisy TAC for background, liver, and tumor.We have performed two sets of experiments. For these three experiments, we set an error margin forCpn and compare the propagated errors in k1, k2, k3, and k4 of liver and tumor using the derived expressions and optimization. Optimization is a numerical method using the conjugate gradient method (CGM).First, we test the mathematical model in one dimension and assume that one of the 19 samples ofCpn has an error −20% to 20% defined as (Cpn(error)-Cpn(true))/Cpn(true). Figure 4 shows a comparison of the results from the derived expressions and optimization for two arbitrary samples of Cp. We observe that the derived expressions provide very accurate approximations of the the errors in k1, k2, k3, and k4 of liver.Comparison between the estimatedk1, k2, k3, and k4 of liver using the derived expressions and numerical method for a range of erroneous blood functions. A single sample out of 19 samples of Cpn has an error. The results are given for two random samples.
(a)(b)(c)(d)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. Figure5 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.Figure 5
Three compartments used to model the transfer of the tracer between physical compartments and chemical states.Figures6 and 7 show that the results from the derived expressions are very close to ones obtained from numerical optimization.Comparison between the estimatedk1, k2, k3, and k4 of liver using the derived expressions and numerical approximation for a range of erroneous blood functions. Top two figures show results when all samples of the blood are erroneous, (a,b)/(c,d) when the initial peak is erroneous, and (e,f) when the tail part is erroneous.
(a)(b)(c)(d)(e)(f)Comparison between the estimatedk1, k2, k3, and k4 of tumor using the derived expressions and numerical approximation for a range of erroneous blood functions. Top two figures show results when all samples of the blood are erroneous, middle two when the initial peak is erroneous, and the bottom two when the tail part is erroneous.
(a)(b)(c)(d)(e)(f)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 ink2, k3, and k4 are negligibly small while change in k1 is relatively large. k1 of liver drops from 1.0093 to 0.6728, k1 of background drops from 0.0126 to 0.0084, k1 of tumor drops from 0.7599 to 0.5067. And k2, k3, and k4 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 thatk1 and k3 deviate from the true value considerably, k4 also deviates from the true value but not as much, while k2 almost remains the same, indicating that the error in the initial peak affects the estimation of kinetic parameter k1, k3, and k4 more than k2. k1 of liver decreases from 1.0122 to 0.6766, k2 of liver increases from 0.9784 to 0.9866, k3 of liver increases from 0.0110 to 0.0129, k4 of liver decreases from 0.0195 to 0.0155. k1 of background decreases from 0.0132 to 0.0081, k2 of background decreases from 0.0594 to 0.0433, k3 of background increases from 0.0001 to 0.0017, k4 of background increases from −0.0215 to 0.0055. k1 of tumor decreases from 0.7455 to 0.5142, k2 of tumor increases from 0.6545 to 0.7124, k3 of tumor increases from 0.0449 to 0.0552, and k4 of tumor decreases from 0.0024 to 0.0003.(c)
For Case III, we observe a different effect; the error in the tail part affects the estimation of kinetic parameterk3 most, k4 to a lesser degree but does not affect k1 and k2. k1 of liver decreases from 0.8141 to 0.8118, k2 of liver decreases from 0.9941 to 0.9815, k3 of liver decreases from 0.0156 to 0.0095, k4 of liver decreases from 0.0186 to 0.0155. k1 of background increases from 0.0097 to 0.0105, k2 of background increases from 0.0397 to 0.0601, k3 of background increases from 0.0008922 to 0.0009358, k4 of background decreases from 0.0025 to −0.0013. k1 of tumor decreases from 0.6176 to 0.6029, k2 of tumor decreases from 0.7173 to 0.6690, k3 of tumor decreases from 0.0588 to 0.0446, and k4 of tumor increases from 0.0005 to 0.0018.In addition to the figures, Table1 lists the relations between the error rate of Cp and the percentage of changes in estimation of the parameters k1, k2, k3, and k4 of the liver in three cases defined before.Table 1
The percentage error in estimatedk1, k2, k3, and k4 of liver affected by a range of erroneous blood functions in three cases.
Error rate ofCpPercentage error in parametersCase Ik1k2k3k410%10%~0~0~020%20%~0~0~0Case IIk1k2k3k410%10%~012.5%4%20%20%~025%8%Case IIIk1k2k3k410%~0~010%5%20%~0~020%10%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 parameterk1, k3, and k4 more than k2. The parameter k3 changes 12.5% when error rate of Cp becomes 10% while k1 changes 10% and k4 changes 4%.(b)
And the error in the tail part affects the estimation of kinetic parameterk3 most, k4 second but does not affect k1 and k2. The parameter k3 will change 10% when error rate of Cp becomes 10% while k1 will change 5%.(c)
If the overall blood input function has almost the same error rate, we find that the error in the parametersk2, k3, and k4 is negligibly small while k1 changes relatively large. In this case, only the parameter k1 will change 10% when error rate of Cp becomes 10%.Table2  lists the relations between the error rate of Cp and the percentage of changes in estimation of the parameters k1, k2, k3, and k4 of the tumor in three cases defined before. Similar conclusion can be drawn from the data of the table.Table 2
The percentage error in estimatedk1, k2, k3, and k4 of tumor affected by a range of erroneous blood functions in three cases.
Error rate of CpPercentage error in parametersCase Ik1k2k3k410%10%~0~0~020%20%~0~0~0Case IIk1k2k3k410%10%~07.5%25%20%20%~015%50%Case IIIk1k2k3k410%~0~09%30%20%~0~018%60%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.


## 4. 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 ink1, k2, k3, and k4 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 k1, k3, and k4 more than k2. The error in the tail part affects the estimation of kinetic parameter k3 most, k4 second but does not affect k1 and k2. If the overall blood input function has almost the same error rate, we find that the error in k2, k3, and k4 is negligibly small while k1 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.



---

*Source: 234679-2010-11-28.xml*