Modern Mathematics in Biomedical ImagingView this Special Issue
Research Article | Open Access
Yafang Cheng, İmam Şamil Yetik, "Propagation of Blood Function Errors to the Estimates of Kinetic Parameters with Dynamic PET", International Journal of Biomedical Imaging, vol. 2011, Article ID 234679, 12 pages, 2011. https://doi.org/10.1155/2011/234679
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.
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 .
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 ) and photodynamic therapy . 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 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 denote the TAC of a specific tissue, then the relation between and the impulse response of a region can be obtained using diffusion equations.
The differential equations that describe the FDG three-compartment model (Figure 1) are expressed as follows: where denotes the tracer concentration in the plasma assumed to be spatially constant, denotes the concentration of unbounded tracer, is the concentration of the bounded tracer, and denotes the time coordinate. The solution is  for , where “” denotes convolution and We use the following transformations as commonly done  since the transformed parameters are more suitable for our mathematical model: with inverse transforms Then, for the three-compartment tissue modeling, the relation between and the impulse response containing the kinetic parameters is where is
Estimation of the kinetic parameters , , , and for three-compartment tissue modeling based on requires that the blood function 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 , , , and 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 times for noisy TAC measurements we can use the following model which can be written as where is the convolution matrix of the impulse response of the tissue for region , denotes the vector of the blood function, and 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 . 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  and the Runge-Kutte methods .
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 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.
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 parameters , , , are Our goal is to obtain the errors in the kinetic parameters , , , and stemming from the error in the blood function . The first step in this calculation is to calculate the derivatives of the implicit estimator with respect to the elements of at a fixed 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 .
For a solution 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 at a fixed point that maps the into an estimate . By applying the chain rule to this equation, we can differentiate the above two equations with respect to and obtain and similar equations for , , and can be obtained easily.
Thus we derive the derivative expression that we are interested in, Because of the presence of , ,, and in the terms , , , and , we cannot simply integrate , , , and to obtain the errors in the kinetic parameter estimates. However, we can calculate , , , and at a fixed point of , and a sequential procedure can be applied to calculate a new value of , , , and , and we can use them to evaluate , , , and at the next value until the complete range of is covered. This procedure can be performed by methods such as Runge-Kutte methods, predictor-corrector method and Richardson extrapolation. where , for , , , and , and is the small step size we define according to the demand on accuracy and speed.
This completes the calculation of the errors in , , , and stemming from the error in estimation of . To summarize the expressions derived in this section provide the errors in the kinetic parameters , , , stemming from errors in the blood . 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 , , , and to obtain , , , and . 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 with 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: , , , (, , , ); the second set of background: , , , (, , , ); the third set of tumor: , , , (, , , ). These values are selected to produce a realistic case in . The blood function is produced using the model : with typical values from  The input blood function can be seen in Figure 2.
We simulate a noisy TAC value as follows: where 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 and compare the propagated errors in , , , and 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 of has an error −20% to 20% defined as . Figure 4 shows a comparison of the results from the derived expressions and optimization for two arbitrary samples of . We observe that the derived expressions provide very accurate approximations of the the errors in , , , and 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.
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 , , and are negligibly small while change in is relatively large. of liver drops from 1.0093 to 0.6728, of background drops from 0.0126 to 0.0084, of tumor drops from 0.7599 to 0.5067. And , , and 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 and but inversely scale and . (b)For Case II, we observe that and deviate from the true value considerably, also deviates from the true value but not as much, while almost remains the same, indicating that the error in the initial peak affects the estimation of kinetic parameter , , and more than . of liver decreases from 1.0122 to 0.6766, of liver increases from 0.9784 to 0.9866, of liver increases from 0.0110 to 0.0129, of liver decreases from 0.0195 to 0.0155. of background decreases from 0.0132 to 0.0081, of background decreases from 0.0594 to 0.0433, of background increases from 0.0001 to 0.0017, of background increases from −0.0215 to 0.0055. of tumor decreases from 0.7455 to 0.5142, of tumor increases from 0.6545 to 0.7124, of tumor increases from 0.0449 to 0.0552, and 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 parameter most, to a lesser degree but does not affect and . of liver decreases from 0.8141 to 0.8118, of liver decreases from 0.9941 to 0.9815, of liver decreases from 0.0156 to 0.0095, of liver decreases from 0.0186 to 0.0155. of background increases from 0.0097 to 0.0105, of background increases from 0.0397 to 0.0601, of background increases from 0.0008922 to 0.0009358, of background decreases from 0.0025 to −0.0013. of tumor decreases from 0.6176 to 0.6029, of tumor decreases from 0.7173 to 0.6690, of tumor decreases from 0.0588 to 0.0446, and of tumor increases from 0.0005 to 0.0018.
In addition to the figures, Table 1 lists the relations between the error rate of and the percentage of changes in estimation of the parameters , , , and 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 , , and more than . The parameter changes 12.5% when error rate of becomes 10% while changes 10% and changes 4%. (b)And the error in the tail part affects the estimation of kinetic parameter most, second but does not affect and . The parameter will change 10% when error rate of becomes 10% while will change 5%. (c)If the overall blood input function has almost the same error rate, we find that the error in the parameters , , and is negligibly small while changes relatively large. In this case, only the parameter will change 10% when error rate of becomes 10%.
Table 2 lists the relations between the error rate of and the percentage of changes in estimation of the parameters , , , and 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.
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 , , , and 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 , , and more than . The error in the tail part affects the estimation of kinetic parameter most, second but does not affect and . If the overall blood input function has almost the same error rate, we find that the error in , , and is negligibly small while 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.
For a solution 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 that maps the into an estimate . We can rewrite (13) as By applying the chain rule to this equation, we can differentiate the above two equations with respect to and obtain and similar equations for , , and 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 where “” denotes the transpose. The terms , ,, , , , , and can be calculated using the compartment model where “” denotes an operation converting a vector to its corresponding convolution matrix.
In this paper, we modify regular Runge-Kutte methods  for ODE to adapt to partial differential equations. By far the most common approximation is the fourth-order Runge-Kutte approximation. where , for , , , and , is the small step size we define according to the demand on accuracy and speed, and is a scalar given by the determinant
Let us define where is as defined before. In multidimension, we perform the calculations times sequentially for every for a single step
- M. N. Wernick and J. N. Aarsvold, Eds., Emission Tomography: The Fundamentals of PET and SPECT, Elsevier, New York, NY, USA, 2004.
- J. X. Ji and Y. Jiraraksopakun, “Model-based simulation of dynamic magnetic resonance imaging signals,” Biomedical Signal Processing and Control, vol. 3, no. 4, pp. 305–311, 2008.
- T. Bastogne, L. Tirand, D. Bechet, M. Barberi-Heyob, and A. Richard, “System identification of photosensitiser uptake kinetics in photodynamic therapy,” Biomedical Signal Processing and Control, vol. 2, no. 3, pp. 217–225, 2007.
- M. E. Phelps, S. C. Huang, and E. J. Hoffman, “Tomographic measurement of local cerebral glucose metabolic rate in humans with (F-18)2-fluoro-2-deoxy-D-glucose: validation of method,” Annals of Neurology, vol. 6, no. 5, pp. 371–388, 1979.
- A. A. Lammertsma and S. P. Hume, “Simplified reference tissue model for PET receptor studies,” NeuroImage, vol. 4, no. 3, pp. 153–158, 1996.
- H. Watabe, R. E. Carson, and H. Iida, “The reference tissue model : three compartments for the reference region,” NeuroImage, vol. 11, no. 6, p. S12, 2000.
- J. Tseng, L. K. Dunnwald, E. K. Schubert et al., “18F-FDG kinetics in locally advanced breast cancer: correlation with tumor blood flow and changes in response to neoadjuvant chemotherapy,” Journal of Nuclear Medicine, vol. 45, no. 11, pp. 1829–1837, 2004.
- S. Chen, C. Ho, D. Feng, and Z. Chi, “Tracer kinetic modeling of 11C-acetate applied in the liver with positron emission tomography,” IEEE Transactions on Medical Imaging, vol. 23, no. 4, pp. 426–432, 2004.
- I. Ş. Yetik and J. Y. Qi, “Direct estimation of kinetic parameters from the sinogram with an unknown blood function,” in Proceedings of the 3rd IEEE International Symposium on Biomedical Imaging, pp. 295–298, April 2006.
- D. Y. Riabkov and E. V. R. Di Bella, “Estimation of kinetic parameters without input functions: analysis of three methods for multichannel blind identification,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 11, pp. 1318–1327, 2002.
- M. I. Gürelli and C. L. Nikias, “EVAM: an eigenvector-based algorithm for multichannel blind deconvolution of input colored signals,” IEEE Transactions on Signal Processing, vol. 43, no. 1, pp. 134–149, 1995.
- D. Feng, K.-P. Wong, C.-M. Wu, and W.-C. Siu, “A technique for extracting physiological parameters and the required input fuction smultaneously from PET iage masurements: theory and smulation study,” IEEE Transactions on Information Technology in Biomedicine, vol. 1, no. 4, pp. 243–254, 1997.
- K.-P. Wong, D. Feng, S. R. Meikle, and M. J. Fulham, “Simultaneous estimation of physiological parameters and the input function—In vivo PET data,” IEEE Transactions on Information Technology in Biomedicine, vol. 5, no. 1, pp. 67–76, 2001.
- Y. Cheng and I. Ş. Yetik, “Effect of the blood function error on the estimated kinetic parameters with dynamic PET,” in Proceedings of the 5th IEEE International Symposium on Biomedical Imaging (ISBI '08), pp. 1585–1588, May 2008.
- J. A. Fessier, “Mean and variance of implicitly defined biased estimators (Such as penalized maximum likelihood): applications to tomography,” IEEE Transactions on Image Processing, vol. 5, no. 3, pp. 493–506, 1996.
- N. Gershenfeld, The Nature of Mathematical Modeling, Cambridge University Press, Cambridge, UK, 1999.
Copyright © 2011 Yafang Cheng and İmam Şamil Yetik. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.