Home | About | Journals | Submit | Contact Us | Français |

**|**Int J Biomed Imaging**|**v.2011; 2011**|**PMC2993041

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Derivation of the Error Propagation for Three-Compartment Tissue Modeling
- 3. Computer Simulations
- 4. Conclusion
- References

Authors

Related links

Int J Biomed Imaging. 2011; 2011: 234679.

Published online 2010 November 28. doi: 10.1155/2011/234679

PMCID: PMC2993041

Department of Electrical Engineering, Medical Imaging Research Center, Illinois Institute of Technology, Chicago, IL 60616-3793, USA

*İmam Şamil Yetik: Email: ude.tii@kitey

Academic Editor: Robert J. Plemmons

Received 2010 May 27; Accepted 2010 August 16.

Copyright © 2011 Y. Cheng and İ. Ş. 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.

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 [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 *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 three-compartment model (Figure 1) are expressed as follows:

Three compartments used to model the transfer of the tracer between physical compartments and chemical states.

$$\begin{array}{cc}\hfill \frac{\text{d}{C}_{u}\left(t\right)}{\text{d}t}& ={k}_{1}{C}_{p}\left(t\right)-({k}_{2}+{k}_{3}){C}_{u}\left(t\right)+{k}_{4}{C}_{b}\left(t\right),\hfill \\ \hfill \frac{\text{d}{C}_{b}\left(t\right)}{\text{d}t}& ={k}_{3}{C}_{u}\left(t\right)-{k}_{4}{C}_{b}\left(t\right),\hfill \end{array}$$

(1)

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]

$$\begin{array}{cc}\hfill {C}_{u}\left(t\right)& =\frac{{k}_{1}}{c-d}\left[\left({k}_{4}-d\right){e}^{-dt}+\left(c-{k}_{4}\right){e}^{-ct}\right]\ast {C}_{p}\left(t\right),\hfill \\ \hfill {C}_{b}\left(t\right)& =\frac{{k}_{1}{k}_{3}}{c-d}\left[{e}^{-dt}-{e}^{-ct}\right]\ast {C}_{p}\left(t\right),\hfill \end{array}$$

(2)

for *t* > 0, where “” denotes convolution and

$$\begin{array}{c}d=\frac{\left({k}_{2}+{k}_{3}+{k}_{4}\right)}{2}-\frac{\sqrt{{\left({k}_{2}+{k}_{3}+{k}_{4}\right)}^{2}-4{k}_{2}{k}_{4}}}{2},\\ c=\frac{\left({k}_{2}+{k}_{3}+{k}_{4}\right)}{2}+\frac{\sqrt{{\left({k}_{2}+{k}_{3}+{k}_{4}\right)}^{2}-4{k}_{2}{k}_{4}}}{2}.\end{array}$$

(3)

We use the following transformations as commonly done [9] since the transformed parameters are more suitable for our mathematical model:

$$\begin{array}{cc}\hfill a& =\frac{{k}_{1}}{2\mathrm{\Delta}}\left({k}_{2}-{k}_{3}-{k}_{4}+\mathrm{\Delta}\right),\hfill \\ \hfill b& =\frac{{k}_{1}}{2\mathrm{\Delta}}\left(-{k}_{2}+{k}_{3}+{k}_{4}+\mathrm{\Delta}\right),\hfill \\ \hfill c& =0.5\left({k}_{2}+{k}_{3}+{k}_{4}+\mathrm{\Delta}\right),\hfill \\ \hfill d& =0.5\left({k}_{2}+{k}_{3}+{k}_{4}-\mathrm{\Delta}\right),\hfill \\ \hfill \mathrm{\Delta}& =\left|\sqrt{{\left({k}_{2}+{k}_{3}+{k}_{4}\right)}^{2}-4{k}_{2}{k}_{4}}\right|\hfill \end{array}$$

(4)

with inverse transforms

$$\begin{array}{cc}\hfill {k}_{1}& =a+b,\hfill \\ \hfill {k}_{2}& =\frac{ac+bd}{a+b},\hfill \\ \hfill {k}_{3}& =\frac{ab{\left(c-d\right)}^{2}}{\left(a+b\right)\left(ac+bd\right)},\hfill \\ \hfill {k}_{4}& =\frac{cd\left(a+b\right)}{ac+bd}.\hfill \end{array}$$

(5)

Then, for the three-compartment tissue modeling, the relation between *f*(*t*) and the impulse response *h*(*t*) containing the kinetic parameters is

$$\begin{array}{c}f\left(t\right)=h\left(t\right)\ast {C}_{p}\left(t\right),\end{array}$$

(6)

where *h*(*t*) is

$$\begin{array}{c}h\left(t\right)=a{e}^{-ct}+b{e}^{-dt}.\end{array}$$

(7)

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 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–13] and references therein.

Since we have discrete measurements at times *t* = Δ*t* · *n* for noisy TAC measurements we can use the following model

$$\begin{array}{c}{f}^{\left(i\right)}\left(n\right)={h}^{\left(i\right)}\left(n\right)\ast {C}_{p}\left(n\right)+{\epsilon}^{\left(i\right)}\left(n\right),\end{array}$$

(8)

which can be written as

$$\begin{array}{c}{\stackrel{\u20d7}{f}}^{\left(i\right)}={H}^{\left(i\right)}\stackrel{\u20d7}{{C}_{p}}+{\stackrel{\u20d7}{\epsilon}}^{\left(i\right)},\end{array}$$

(9)

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 ${\stackrel{\u20d7}{\epsilon}}^{(i)}$ is the noise vector. Stacking different regions of interest together, we can write the equations in the following form

$$\begin{array}{c}\stackrel{\u20d7}{F}=H\stackrel{\u20d7}{{C}_{p}}+\stackrel{\u20d7}{\epsilon}.\end{array}$$

(10)

We can estimate the kinetic parameters and the blood function by minimizing the following cost function:

$$\begin{array}{c}R={||\stackrel{\u20d7}{F}-\widehat{H}\stackrel{\u20d7}{{\widehat{C}}_{p}}||}^{2}.\end{array}$$

(11)

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 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.

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 $\widehat{a}$, $\widehat{b}$, $\widehat{c}$, $\widehat{d}$ are

$$\begin{array}{c}\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right]=\mathrm{arg}\hspace{0.17em}\underset{\widehat{a},\widehat{b},\widehat{c},\widehat{d}}{\mathrm{min}\hspace{0.17em}\hspace{0.17em}}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right).\end{array}$$

(12)

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

$${0=\frac{\partial}{\partial a}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right)|}_{a=\widehat{a}}$$

(13)

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

$$\begin{array}{c}\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right]=g\left({\widehat{C}}_{p}\right)={\left[{g}_{1}\left({\widehat{C}}_{p}\right),{g}_{2}\left({\widehat{C}}_{p}\right),{g}_{3}\left({\widehat{C}}_{p}\right),{g}_{4}\left({\widehat{C}}_{p}\right)\right]}^{\prime}\end{array}$$

(14)

that maps the ${\widehat{C}}_{p}$ into an estimate $[\widehat{a},\widehat{b},\widehat{c},\widehat{d}]$. By applying the chain rule to this equation, we can differentiate the above two equations with respect to ${\widehat{C}}_{p}$ and obtain

$$\begin{array}{cc}\hfill 0& =\frac{{\partial}^{2}}{\partial {a}^{2}}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{1}\left({\widehat{C}}_{p}\right)\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial b}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{2}\left({\widehat{C}}_{p}\right)\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial c}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{3}\left({\widehat{C}}_{p}\right)\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial d}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{4}\left({\widehat{C}}_{p}\right)\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{p}}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right),\hfill \end{array}$$

(15)

and similar equations for *b*, *c*, and *d* can be obtained easily.

Thus we derive the derivative expression that we are interested in,

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\frac{\partial {g}_{1}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \\ \frac{\partial {g}_{2}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \\ \frac{\partial {g}_{3}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \\ \frac{\partial {g}_{4}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\end{array}\right]& ={\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right]}^{-1}\hfill \\ \hfill & \hspace{1em}\times \left[\begin{array}{c}-\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R\\ \\ -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R\\ \\ -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R\\ \\ -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R\end{array}\right].\hfill \end{array}$$

(16)

Because of the presence of , *b*,*c*, and *d* in the terms $\partial \text{a}/\partial {\widehat{C}}_{{p}_{n}}$, $\partial b/\partial {\widehat{C}}_{{p}_{n}}$, $\partial c/\partial {\widehat{C}}_{{p}_{n}}$, and $\partial d/\partial {\widehat{C}}_{{p}_{n}}$, we cannot simply integrate $(\partial \text{a}/\partial {\widehat{C}}_{{p}_{n}})\Delta {\widehat{C}}_{{p}_{n}}$, $(\partial b/\partial {\widehat{C}}_{{p}_{n}})\Delta {\widehat{C}}_{{p}_{n}}$, $(\partial c/\partial {\widehat{C}}_{{p}_{n}})\Delta {\widehat{C}}_{{p}_{n}}$, and $(\partial d/\partial {\widehat{C}}_{{p}_{n}})\Delta {\widehat{C}}_{{p}_{n}}$ to obtain the errors in the kinetic parameter estimates. However, we can calculate $\partial \text{a}/\partial {\widehat{C}}_{{p}_{n}}$, $\partial b/\partial {\widehat{C}}_{{p}_{n}}$, $\partial c/\partial {\widehat{C}}_{{p}_{n}}$, and $\partial d/\partial {\widehat{C}}_{{p}_{n}}$ at a fixed point of ${\widehat{C}}_{{p}_{n}}$, 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 $\partial \text{a}/\partial {\widehat{C}}_{{p}_{n}}$, $\partial b/\partial {\widehat{C}}_{{p}_{n}}$, $\partial c/\partial {\widehat{C}}_{{p}_{n}}$, and $\partial d/\partial {\widehat{C}}_{{p}_{n}}$ at the next *C*_{p} value until 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.

$$\begin{array}{cc}\hfill {s}_{1,i}& =h{u}_{i}\left({C}_{{p}_{n}},a,b,c,d\right),\hfill \\ \hfill {s}_{2,i}& =h{u}_{i}\left({C}_{{p}_{n}}+\frac{h}{2},a+\frac{{s}_{\mathrm{1,1}}}{2},b+\frac{{s}_{\mathrm{1,2}}}{2},c+\frac{{s}_{\mathrm{1,3}}}{2},d+\frac{{s}_{\mathrm{1,4}}}{2}\right),\hfill \\ \hfill {s}_{3,i}& =h{u}_{i}\left({C}_{{p}_{n}}+\frac{h}{2},a+\frac{{s}_{\mathrm{2,1}}}{2},b+\frac{{s}_{\mathrm{2,2}}}{2},c+\frac{{s}_{\mathrm{2,3}}}{2},d+\frac{{s}_{\mathrm{2,4}}}{2}\right),\hfill \\ \hfill {s}_{4,i}& =h{u}_{i}\left({C}_{{p}_{n}}+h,a+{s}_{\mathrm{3,1}},b+{s}_{\mathrm{3,2}},c+{s}_{\mathrm{3,3}},d+{s}_{\mathrm{3,4}}\right),\hfill \\ \hfill {y}_{i}\left({C}_{{p}_{n}}+h\right)& ={y}_{i}\left({C}_{{p}_{n}}\right)+\frac{{s}_{1,i}}{6}+\frac{{s}_{2,i}}{3}+\frac{{s}_{3,i}}{3}+\frac{{s}_{4,\text{i}}}{6},\hfill \end{array}$$

(17)

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.

We apply the proposed mathematical model to simulated dynamic PET data with *t* = [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: *k*_{1} = 0.8, *k*_{2} = 0.98, *k*_{3} = 0.012, *k*_{4} = 0.017 (*a* = 0.79, *b* = 0.01, *c* = 0.9922, *d* = 0.0168); the second set of background: *k*_{1} = 0.01, *k*_{2} = 0.05, *k*_{3} = 0.001, *k*_{4} = 0.0001 (*a* = 0.0098, *b* = 0.00019683, *c* = 0.0510, *d* = 0.000098035); the third set of tumor: *k*_{1} = 0.598, *k*_{2} = 0.680, *k*_{3} = 0.05, *k*_{4} = 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]:

$$\begin{array}{c}{C}_{p}\left(t\right)=({b}_{1}t-{b}_{2}-{b}_{3}){e}^{{\lambda}_{1}t}+{b}_{2}{e}^{{\lambda}_{2}t}+{b}_{3}{e}^{{\lambda}_{3}t}\end{array}$$

(18)

with typical values from [12]

$$\left[{b}_{1},{b}_{2},{b}_{3},{\lambda}_{1},{\lambda}_{2},{\lambda}_{3}\right]=\left[\mathrm{900,22,21},-4.1,-0.12,-0.01\right]{\mathrm{min}\hspace{0.17em}}^{-1}.$$

(19)

The input blood function can be seen in Figure 2.

We simulate a noisy TAC value as follows:

$$\begin{array}{c}{F}_{\text{noisy}}={F}_{\text{original}}\xb7(1+N),\end{array}$$

(20)

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).

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.

Comparison between the estimated *k*_{1}, *k*_{2}, *k*_{3}, and *k*_{4} of liver using the derived expressions and numerical method for a range of erroneous blood functions. A single sample out of 19 samples of *C*_{pn} has an error. The results are given for two random samples. **...**

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.

Three compartments used to model the transfer of the tracer between physical compartments and chemical states.

Figures Figures66 and and77 show that the results from the derived expressions are very close to ones obtained from numerical optimization.

Comparison between the estimated *k*_{1}, *k*_{2}, *k*_{3}, and *k*_{4} 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 **...**

Comparison between the estimated *k*_{1}, *k*_{2}, *k*_{3}, and *k*_{4} 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 following items summarize the effect of the blood function error to the final kinetic parameter estimates.

- 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 Figures5,5, ,6,6, and and77 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*. - 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}.*k*_{1}of liver decreases from 1.0122 to 0.6766,*k*_{2}of liver increases from 0.9784 to 0.9866,*k*_{3}of liver increases from 0.0110 to 0.0129,*k*_{4}of liver decreases from 0.0195 to 0.0155.*k*_{1}of background decreases from 0.0132 to 0.0081,*k*_{2}of background decreases from 0.0594 to 0.0433,*k*_{3}of background increases from 0.0001 to 0.0017,*k*_{4}of background increases from −0.0215 to 0.0055.*k*_{1}of tumor decreases from 0.7455 to 0.5142,*k*_{2}of tumor increases from 0.6545 to 0.7124,*k*_{3}of tumor increases from 0.0449 to 0.0552, and*k*_{4}of tumor decreases from 0.0024 to 0.0003. - For Case III, we observe a different effect; the error in the tail part affects the estimation of kinetic parameter
*k*_{3}most,*k*_{4}to a lesser degree but does not affect*k*_{1}and*k*_{2}.*k*_{1}of liver decreases from 0.8141 to 0.8118,*k*_{2}of liver decreases from 0.9941 to 0.9815,*k*_{3}of liver decreases from 0.0156 to 0.0095,*k*_{4}of liver decreases from 0.0186 to 0.0155.*k*_{1}of background increases from 0.0097 to 0.0105,*k*_{2}of background increases from 0.0397 to 0.0601,*k*_{3}of background increases from 0.0008922 to 0.0009358,*k*_{4}of background decreases from 0.0025 to −0.0013.*k*_{1}of tumor decreases from 0.6176 to 0.6029,*k*_{2}of tumor decreases from 0.7173 to 0.6690,*k*_{3}of tumor decreases from 0.0588 to 0.0446, and*k*_{4}of tumor increases from 0.0005 to 0.0018.

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.

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.

Our conclusion of the relation between the blood function error and the error on the kinetic parameters can be summarized as follows:

- 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%. - 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%. - 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.

The percentage error in estimated *k*_{1}, *k*_{2}, *k*_{3}, and *k*_{4} of tumor affected by a range of erroneous blood functions in three cases.

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 *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.

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

$$\begin{array}{c}{0=\frac{\partial}{\partial a}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right)|}_{a=\widehat{a}},\\ {0=\frac{\partial}{\partial b}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right)|}_{b=\widehat{b}},\\ {0=\frac{\partial}{\partial c}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right)|}_{c=\widehat{c}},\\ {0=\frac{\partial}{\partial d}R\left(\left[\widehat{a},\widehat{b},\widehat{c},\widehat{d}\right],{\widehat{C}}_{p}\right)|}_{d=\widehat{d}}.\end{array}$$

(A.1)

Let us define an implicit function which will be convenient for the application of chain rule

$$\begin{array}{c}[\widehat{a},\widehat{b},\widehat{c},\widehat{d}]=g\left({\widehat{C}}_{p}\right)={\left[{g}_{1}\left({\widehat{C}}_{p}\right),{g}_{2}\left({\widehat{C}}_{p}\right),{g}_{3}\left({\widehat{C}}_{p}\right),{g}_{4}\left({\widehat{C}}_{p}\right)\right]}^{\prime}\end{array}$$

(A.2)

that maps the ${\widehat{C}}_{p}$ into an estimate $[\widehat{a},\widehat{b},\widehat{c},\widehat{d}]$. We can rewrite (13) as

$$\begin{array}{c}0=\frac{\partial}{\partial a}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right),\\ 0=\frac{\partial}{\partial b}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right),\\ 0=\frac{\partial}{\partial c}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right),\\ 0=\frac{\partial}{\partial d}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right).\end{array}$$

(A.3)

By applying the chain rule to this equation, we can differentiate the above two equations with respect to ${\widehat{C}}_{p}$ and obtain

$$\begin{array}{cc}\hfill 0& =\frac{{\partial}^{2}}{\partial {a}^{2}}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{1}\left({\widehat{C}}_{p}\right),\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial b}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{2}\left({\widehat{C}}_{p}\right),\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial c}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{3}\left({\widehat{C}}_{p}\right),\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial d}R\left(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}\right)\frac{\partial}{\partial {\widehat{C}}_{p}}{g}_{4}\left({\widehat{C}}_{p}\right),\hfill \\ \hfill & \hspace{1em}+\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{p}}R(g\left({\widehat{C}}_{p}\right),{\widehat{C}}_{p}),\hfill \end{array}$$

(A.4)

and similar equations for *b*, *c*, and *d* can be obtained easily.

We can write these four equations in matrix form

$$\begin{array}{cc}\hfill \left[\begin{array}{c}-\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R\end{array}\right]& =\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right]\hfill \\ \hfill & \hspace{1em}\times \left[\begin{array}{c}\frac{\partial {g}_{1}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{2}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{3}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{4}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\end{array}\right].\hfill \end{array}$$

(A.5)

Then a simple matrix inversion provides

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\frac{\partial {g}_{1}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{2}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{3}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\\ \frac{\partial {g}_{4}\left({\widehat{C}}_{p}\right)}{\partial {\widehat{C}}_{{p}_{n}}}\end{array}\right]& ={\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right]}^{-1}\hfill \\ \hfill & \hspace{1em}\times \left[\begin{array}{c}-\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R\\ -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R\end{array}\right].\hfill \end{array}$$

(A.6)

The elements of the matrix and the vector on the right hand side can be calculated based on the cost function

$$\begin{array}{cc}\hfill \frac{{\partial}^{2}R}{\partial {a}^{2}}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial a}\right)}^{\prime}\frac{\partial H}{\partial a}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial a\partial b}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial a}\right)}^{\prime}\frac{\partial H}{\partial b}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial a\partial c}& =-2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial a\partial c}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial a\partial c}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial a}\right)}^{\prime}\frac{\partial H}{\partial c}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial a\partial d}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial a}\right)}^{\prime}\frac{\partial H}{\partial d}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial {b}^{2}}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial b}\right)}^{\prime}\frac{\partial H}{\partial b}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial b\partial c}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial b}\right)}^{\prime}\frac{\partial H}{\partial c}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial b\partial d}& =-2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial b\partial d}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial b\partial d}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}\left(\frac{\partial H}{\partial b}\right)\frac{\partial H}{\partial d}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial {c}^{2}}& =-2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial {c}^{2}}\right)}^{\prime}\stackrel{\u20d7}{F}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial {c}^{2}}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial c}\right)}^{\prime}\frac{\partial H}{\partial c}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial c\partial d}& =2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial c}\right)}^{\prime}\frac{\partial H}{\partial d}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial {d}^{2}}& =-2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial {d}^{2}}\right)}^{\prime}\stackrel{\u20d7}{F}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{{\partial}^{2}H}{\partial {d}^{2}}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial d}\right)}^{\prime}\frac{\partial H}{\partial d}\stackrel{\u20d7}{{C}_{p}},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial a\partial {\widehat{C}}_{{p}_{n}}}& =-2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial a}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial a}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial a}\right)}^{\prime}H{\left[0\cdots 0{1}_{i=n}0\cdots 0\right]}^{\prime},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial b\partial {\widehat{C}}_{{p}_{n}}}& =-2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial b}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial b}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial b}\right)}^{\prime}H{\left[0\cdots 0{1}_{i=n}0\cdots 0\right]}^{\prime},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial c\partial {\widehat{C}}_{{p}_{n}}}& =-2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial c}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial c}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial c}\right)}^{\prime}H{\left[0\cdots 0{1}_{i=n}0\cdots 0\right]}^{\prime},\hfill \\ \hfill \frac{{\partial}^{2}R}{\partial d\partial {\widehat{C}}_{{p}_{n}}}& =-2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial d}\right)}^{\prime}\stackrel{\u20d7}{F}\hfill \\ \hfill & \hspace{1em}+2\left[0\cdots 0{1}_{i=n}0\cdots 0\right]{\left(\frac{\partial H}{\partial d}\right)}^{\prime}H\stackrel{\u20d7}{{C}_{p}}\hfill \\ \hfill & \hspace{1em}+2{\stackrel{\u20d7}{{C}_{p}}}^{\prime}{\left(\frac{\partial H}{\partial d}\right)}^{\prime}H{\left[0\cdots 0{1}_{i=n}0\cdots 0\right]}^{\prime},\hfill \end{array}$$

(A.7)

where “′” denotes the transpose. The terms *H*/*a*, *H*/*b*,*H*/*c*, *H*/*d*, (^{2}/*a**c*)*H*, (^{2}/*c*^{2})*H*, (^{2}/*b**d*)*H*, and (^{2}/*d*^{2})*H* can be calculated using the compartment model

$$\begin{array}{cc}\hfill \frac{\partial H}{\partial a}& =\text{T}\left[{e}^{-ct}\right],\hfill \\ \hfill \frac{\partial H}{\partial b}& =\text{T}\left[{e}^{-dt}\right],\hfill \\ \hfill \frac{\partial H}{\partial c}& =\text{T}[-at{e}^{-ct}],\hfill \\ \hfill \frac{\partial H}{\partial d}& =\text{T}[-bt{e}^{-dt}],\hfill \\ \hfill \frac{{\partial}^{2}}{\partial a\partial c}H& =\text{T}\left[-t{e}^{-ct}\right],\hfill \\ \hfill \frac{{\partial}^{2}}{\partial b\partial d}H& =\text{T}[-t{e}^{-dt}],\hfill \\ \hfill \frac{{\partial}^{2}}{\partial {c}^{2}}H& =\text{T}\left[a{t}^{2}{e}^{-ct}\right],\hfill \\ \hfill \frac{{\partial}^{2}}{\partial {d}^{2}}H& =\text{T}\left[b{t}^{2}{e}^{-dt}\right],\hfill \end{array}$$

(A.8)

where “T” denotes an operation converting a vector to its corresponding convolution matrix.

In this paper, we modify regular Runge-Kutte methods [16] for ODE to adapt to partial differential equations. By far the most common approximation is the fourth-order Runge-Kutte approximation.

$$\begin{array}{cc}\hfill {s}_{1,i}& =h{u}_{i}\left({C}_{{p}_{n}},a,b,c,d\right),\hfill \\ \hfill {s}_{2,i}& =h{u}_{i}\left({C}_{{p}_{n}}+\frac{h}{2},a+\frac{{s}_{\mathrm{1,1}}}{2},b+\frac{{s}_{\mathrm{1,2}}}{2},c+\frac{{s}_{\mathrm{1,3}}}{2},d+\frac{{s}_{\mathrm{1,4}}}{2}\right),\hfill \\ \hfill {s}_{3,i}& =h{u}_{i}\left({C}_{{p}_{n}}+\frac{h}{2},a+\frac{{s}_{\mathrm{2,1}}}{2},b+\frac{{s}_{\mathrm{2,2}}}{2},c+\frac{{s}_{\mathrm{2,3}}}{2},d+\frac{{s}_{\mathrm{2,4}}}{2}\right),\hfill \\ \hfill {s}_{4,i}& =h{u}_{i}\left({C}_{{p}_{n}}+h,a+{s}_{\mathrm{3,1}},b+{s}_{\mathrm{3,2}},c+{s}_{\mathrm{3,3}},d+{s}_{\mathrm{3,4}}\right),\hfill \\ \hfill {y}_{i}({C}_{{p}_{n}}+h)& ={y}_{i}\left({C}_{{p}_{n}}\right)+\frac{{s}_{1,i}}{6}+\frac{{s}_{2,i}}{3}+\frac{{s}_{3,i}}{3}+\frac{{s}_{4,i}}{6},\hfill \end{array}$$

(A.9)

where *k*_{i} = *y*_{i}(*C*_{pn}), *i* = 1,2, 3,4 for *a*, *b*, *c*, and *d*, *h* is the small step size we define according to the demand on accuracy and speed,

$$\begin{array}{cc}\hfill {u}_{1}& =\frac{1}{\varphi}\left[\begin{array}{cccc}-\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ & & & \\ -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ & & & \\ -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ & & & \\ -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right],\hfill \\ \hfill {u}_{2}& =\frac{1}{\varphi}\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& -\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial b\partial a}R& -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial c\partial a}R& -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial d\partial a}R& -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right],\hfill \\ \hfill {u}_{3}& =\frac{1}{\varphi}\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& -\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right],\hfill \\ \hfill {u}_{4}& =\frac{1}{\varphi}\left[\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& -\frac{{\partial}^{2}}{\partial a\partial {\widehat{C}}_{{p}_{n}}}R\\ & & & \\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& -\frac{{\partial}^{2}}{\partial b\partial {\widehat{C}}_{{p}_{n}}}R\\ & & & \\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& -\frac{{\partial}^{2}}{\partial c\partial {\widehat{C}}_{{p}_{n}}}R\\ & & & \\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& -\frac{{\partial}^{2}}{\partial d\partial {\widehat{C}}_{{p}_{n}}}R\end{array}\right],\hfill \end{array}$$

(A.10)

and *ϕ* is a scalar given by the determinant

$$\varphi =\left|\begin{array}{cccc}\frac{{\partial}^{2}}{\partial {a}^{2}}R& \frac{{\partial}^{2}}{\partial a\partial b}R& \frac{{\partial}^{2}}{\partial a\partial c}R& \frac{{\partial}^{2}}{\partial a\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial b\partial a}R& \frac{{\partial}^{2}}{\partial {b}^{2}}R& \frac{{\partial}^{2}}{\partial b\partial c}R& \frac{{\partial}^{2}}{\partial b\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial c\partial a}R& \frac{{\partial}^{2}}{\partial c\partial b}R& \frac{{\partial}^{2}}{\partial {c}^{2}}R& \frac{{\partial}^{2}}{\partial c\partial d}R\\ & & & \\ \frac{{\partial}^{2}}{\partial d\partial a}R& \frac{{\partial}^{2}}{\partial d\partial b}R& \frac{{\partial}^{2}}{\partial d\partial c}R& \frac{{\partial}^{2}}{\partial {d}^{2}}R\end{array}\right|,$$

(A.11)

Let us define

$$\begin{array}{cc}\hfill {a}_{i}& ={y}_{1}\left({C}_{{p}_{n}}\right),\hfill \\ \hfill {b}_{i}& ={y}_{2}\left({C}_{{p}_{n}}\right),\hfill \\ \hfill {c}_{i}& ={y}_{3}\left({C}_{{p}_{n}}\right),\hfill \\ \hfill {d}_{i}& ={y}_{4}\left({C}_{{p}_{n}}\right),\hfill \\ \hfill {a}_{i+1}& ={y}_{1}({C}_{{p}_{n}}+h),\hfill \\ \hfill {b}_{i+1}& ={y}_{2}({C}_{{p}_{n}}+h),\hfill \\ \hfill {c}_{i+1}& ={y}_{3}({C}_{{p}_{n}}+h),\hfill \\ \hfill {d}_{i+1}& ={y}_{4}\left({C}_{{p}_{n}}+h\right),\hfill \end{array}$$

(A.12)

where *h* is as defined before. In multidimension, we perform the calculations *N* times sequentially for every *C*_{pn} for a single step

$$\begin{array}{cc}\hfill {a}_{i+1}& ={a}_{i}+{\displaystyle \sum}_{n=1}^{N}\frac{\partial a}{\partial {\widehat{C}}_{{p}_{n}}}\mathrm{\Delta}{\widehat{C}}_{{p}_{n}},\hfill \\ \hfill {b}_{i+1}& ={b}_{i}+{\displaystyle \sum}_{n=1}^{N}\frac{\partial b}{\partial {\widehat{C}}_{{p}_{n}}}\mathrm{\Delta}{\widehat{C}}_{{p}_{n}},\hfill \\ \hfill {c}_{i+1}& ={c}_{i}+{\displaystyle \sum}_{n=1}^{N}\frac{\partial c}{\partial {\widehat{C}}_{{p}_{n}}}\mathrm{\Delta}{\widehat{C}}_{{p}_{n}},\hfill \\ \hfill {d}_{i+1}& ={d}_{i}+{\displaystyle \sum}_{n=1}^{N}\frac{\partial d}{\partial {\widehat{C}}_{{p}_{n}}}\mathrm{\Delta}{\widehat{C}}_{{p}_{n}}.\hfill \end{array}$$

(A.13)

1. Wernick MN, Aarsvold JN, editors. *Emission Tomography: The Fundamentals of PET and SPECT*. New York, NY, USA: Elsevier; 2004.

2. Ji JX, Jiraraksopakun Y. Model-based simulation of dynamic magnetic resonance imaging signals. *Biomedical Signal Processing and Control*. 2008;3(4):305–311.

3. Bastogne T, Tirand L, Bechet D, Barberi-Heyob M, Richard A. System identification of photosensitiser uptake kinetics in photodynamic therapy. *Biomedical Signal Processing and Control*. 2007;2(3):217–225.

4. Phelps ME, Huang SC, Hoffman EJ. 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*. 1979;6(5):371–388. [PubMed]

5. Lammertsma AA, Hume SP. Simplified reference tissue model for PET receptor studies. *NeuroImage*. 1996;4(3):153–158. [PubMed]

6. Watabe H, Carson RE, Iida H. The reference tissue model : three compartments for the reference region. *NeuroImage*. 2000;11(6):p. S12.

7. Tseng J, Dunnwald LK, Schubert EK, 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*. 2004;45(11):1829–1837. [PubMed]

8. Chen S, Ho C, Feng D, Chi Z. Tracer kinetic modeling of 11C-acetate applied in the liver with positron emission tomography. *IEEE Transactions on Medical Imaging*. 2004;23(4):426–432. [PubMed]

9. Yetik IŞ, Qi JY. Direct estimation of kinetic parameters from the sinogram with an unknown blood function. In: Proceedings of the 3rd IEEE International Symposium on Biomedical Imaging; April 2006; pp. 295–298.

10. Riabkov DY, Di Bella EVR. Estimation of kinetic parameters without input functions: analysis of three methods for multichannel blind identification. *IEEE Transactions on Biomedical Engineering*. 2002;49(11):1318–1327. [PubMed]

11. Gürelli MI, Nikias CL. EVAM: an eigenvector-based algorithm for multichannel blind deconvolution of input colored signals. *IEEE Transactions on Signal Processing*. 1995;43(1):134–149.

12. Feng D, Wong K-P, Wu C-M, Siu W-C. 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*. 1997;1(4):243–254. [PubMed]

13. Wong K-P, Feng D, Meikle SR, Fulham MJ. Simultaneous estimation of physiological parameters and the input function—In vivo PET data. *IEEE Transactions on Information Technology in Biomedicine*. 2001;5(1):67–76. [PubMed]

14. Cheng Y, Yetik IŞ. 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); May 2008; pp. 1585–1588.

15. Fessier JA. Mean and variance of implicitly defined biased estimators (Such as penalized maximum likelihood): applications to tomography. *IEEE Transactions on Image Processing*. 1996;5(3):493–506. [PubMed]

16. Gershenfeld N. *The Nature of Mathematical Modeling*. Cambridge, UK: Cambridge University Press; 1999.

Articles from International Journal of Biomedical Imaging are provided here courtesy of **Hindawi Publishing Corporation**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |