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

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.

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 C_𝑝(𝑡) 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: d𝐶_𝑢(𝑡)d𝑡=𝑘₁𝐶_𝑝(𝑡)−𝑘₂+𝑘₃𝐶_𝑢(𝑡)+𝑘₄𝐶_𝑏(𝑡),d𝐶_𝑏(𝑡)d𝑡=𝑘₃𝐶_𝑢(𝑡)−𝑘₄𝐶_𝑏(𝑡), (1) 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 [8] 𝐶_𝑢(𝑡)=𝑘₁𝑘𝑐−𝑑⁴−𝑑𝑒^−𝑑𝑡+𝑐−𝑘₄𝑒^−𝑐𝑡∗𝐶_𝑝(𝑡),𝐶_𝑏(𝑡)=𝑘₁𝑘₃𝑒𝑐−𝑑^−𝑑𝑡−𝑒^−𝑐𝑡∗𝐶_𝑝(𝑡), (2) for 𝑡>0, where “∗” denotes convolution and 𝑘𝑑=²+𝑘₃+𝑘₄2−𝑘₂+𝑘₃+𝑘₄²−4𝑘₂𝑘₄2,𝑘𝑐=²+𝑘₃+𝑘₄2+𝑘₂+𝑘₃+𝑘₄²−4𝑘₂𝑘₄2. (3) We use the following transformations as commonly done [9] since the transformed parameters are more suitable for our mathematical model: 𝑘𝑎=¹𝑘2Δ²−𝑘₃−𝑘₄,𝑘+Δ𝑏=¹2Δ−𝑘₂+𝑘₃+𝑘₄,𝑘+Δ𝑐=0.5₂+𝑘₃+𝑘₄,𝑘+Δ𝑑=0.5₂+𝑘₃+𝑘₄−Δ,Δ=𝑘₂+𝑘₃+𝑘₄²−4𝑘₂𝑘₄ (4) with inverse transforms 𝑘₁𝑘=𝑎+𝑏,₂=𝑎𝑐+𝑏𝑑,𝑘𝑎+𝑏₃=𝑎𝑏(𝑐−𝑑)²(),𝑘𝑎+𝑏)(𝑎𝑐+𝑏𝑑₄=𝑐𝑑(𝑎+𝑏).𝑎𝑐+𝑏𝑑 (5) Then, for the three-compartment tissue modeling, the relation between 𝑓(𝑡) and the impulse response ℎ(𝑡) containing the kinetic parameters is 𝑓(𝑡)=ℎ(𝑡)∗𝐶_𝑝(𝑡), (6) where ℎ(𝑡) is ℎ(𝑡)=𝑎𝑒^−𝑐𝑡+𝑏𝑒^−𝑑𝑡. (7)

Figure 1 

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

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 𝑓^(𝑖)(𝑛)=ℎ^(𝑖)(𝑛)∗𝐶_𝑝(𝑛)+𝜀^(𝑖)(𝑛), (8) which can be written as 𝑓^(𝑖)=𝐻^(𝑖)𝐶_𝑝+𝜀⃗^(𝑖), (9) 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 𝐹=𝐻𝐶_𝑝+𝜀.⃗ (10)

We can estimate the kinetic parameters and the blood function by minimizing the following cost function: 𝑅=𝐹−𝐻𝐶_𝑝². (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.

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 ̂̂𝑑𝑎𝑎,𝑏,𝑐,=argmin_,𝑏̂_,𝑐_,̂𝑑𝑅𝑎,𝑏̂,𝑐,̂𝑑,𝐶_𝑝. (12) 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 [15].

For a solution 𝑎,𝑏,𝑐,𝑑 that minimizes the cost function the partial derivatives of cost function with respect to the kinetic parameters is zero 𝜕0=𝑅̂̂𝑑,𝐶𝜕𝑎𝑎,𝑏,𝑐,_𝑝𝑎=𝑎 (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 𝐶_𝑝(𝑡) at a fixed 𝐶_𝑝(𝑡) point 𝑎,𝑏̂,𝑐,̂𝑑=𝑔𝐶_𝑝=𝑔₁𝐶_𝑝,𝑔₂𝐶_𝑝,𝑔₃𝐶_𝑝,𝑔₄𝐶_𝑝^′ (14) 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 𝜕0=²𝜕𝑎²𝑅𝑔𝐶_𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₁𝐶_𝑝+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑏^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₂𝐶_𝑝+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑐^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₃𝐶_𝑝+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑑^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₄𝐶_𝑝+𝜕²𝜕𝑎𝜕𝐶_𝑝𝑅𝑔𝐶_𝑝,𝐶_𝑝, (15) and similar equations for 𝑏, 𝑐, and 𝑑 can be obtained easily.

Thus we derive the derivative expression that we are interested in, 𝜕𝑔₁𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₂𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₃𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₄𝐶_𝑝𝜕𝐶_𝑝𝑛=𝜕²𝜕𝑎²𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅𝜕𝜕𝑎𝜕𝑑₂𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅𝜕𝜕𝑏𝜕𝑑₂𝑅𝜕𝜕𝑐𝜕𝑎²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝑅𝜕𝜕𝑑𝜕𝑎²𝑅𝜕𝜕𝑑𝜕𝑏²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅⁻¹×−𝜕²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅. (16) Because of the presence of 𝜕, 𝑏,𝑐, and 𝑑 in the terms 𝜕a/𝜕𝐶_𝑝𝑛, 𝜕𝑏/𝜕𝐶_𝑝𝑛, 𝜕𝑐/𝜕𝐶_𝑝𝑛, and 𝜕𝑑/𝜕𝐶_𝑝𝑛, we cannot simply integrate (𝜕a/𝜕𝐶_𝑝𝑛)Δ𝐶_𝑝𝑛, (𝜕𝑏/𝜕𝐶_𝑝𝑛)Δ𝐶_𝑝𝑛, (𝜕𝑐/𝜕𝐶_𝑝𝑛)Δ𝐶_𝑝𝑛, and (𝜕𝑑/𝜕𝐶_𝑝𝑛)Δ𝐶_𝑝𝑛 to obtain the errors in the kinetic parameter estimates. However, we can calculate 𝜕a/𝜕𝐶_𝑝𝑛, 𝜕𝑏/𝜕𝐶_𝑝𝑛, 𝜕𝑐/𝜕𝐶_𝑝𝑛, 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 𝜕a/𝜕𝐶_𝑝𝑛, 𝜕𝑏/𝜕𝐶_𝑝𝑛, 𝜕𝑐/𝜕𝐶_𝑝𝑛, 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.𝑠_1,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛,𝑠,𝑎,𝑏,𝑐,𝑑_2,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ2𝑠,𝑎+^1,12𝑠,𝑏+^1,22𝑠,𝑐+^1,32𝑠,𝑑+^1,42,𝑠_3,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ2𝑠,𝑎+^2,12𝑠,𝑏+^2,22𝑠,𝑐+^2,32𝑠,𝑑+^2,42,𝑠_4,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ,𝑎+𝑠_3,1,𝑏+𝑠_3,2,𝑐+𝑠_3,3,𝑑+𝑠_3,4,𝑦_𝑖𝐶_𝑝𝑛+ℎ=𝑦_𝑖𝐶_𝑝𝑛+𝑠_1,𝑖6+𝑠_2,𝑖3+𝑠_3,𝑖3+𝑠_4,i6, (17) where 𝑘_𝑖=𝑦_𝑖(𝐶_𝑝𝑛), 𝑖=1,2,3,4 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 𝑡=[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: 𝑘₁=0.8, 𝑘₂=0.98, 𝑘₃=0.012, 𝑘₄=0.017 (𝑎=0.79, 𝑏=0.01, 𝑐=0.9922, 𝑑=0.0168); the second set of background: 𝑘₁=0.01, 𝑘₂=0.05, 𝑘₃=0.001, 𝑘₄=0.0001 (𝑎=0.0098, 𝑏=0.00019683, 𝑐=0.0510, 𝑑=0.000098035); the third set of tumor: 𝑘₁=0.598, 𝑘₂=0.680, 𝑘₃=0.05, 𝑘₄=0.001 (𝑎=0.5569, 𝑏=0.0411, 𝑐=0.7301, 𝑑=0.00093142). These values are selected to produce a realistic case in [9]. The blood function is produced using the model [12]:𝐶_𝑝(𝑡)=𝑏₁𝑡−𝑏₂−𝑏₃𝑒^𝜆1𝑡+𝑏₂𝑒^𝜆2𝑡+𝑏₃𝑒^𝜆3𝑡 (18) with typical values from [12]𝑏₁,𝑏₂,𝑏₃,𝜆₁,𝜆₂,𝜆₃=[900,22,21,−4.1,−0.12,−0.01]min⁻¹. (19) The input blood function can be seen in Figure 2.

Figure 2 

The input blood function.

We simulate a noisy TAC value as follows: 𝐹_noisy=𝐹_original⋅(1+𝑁), (20) 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.

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 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 (𝐶_𝑝𝑛(error)−𝐶_𝑝𝑛(true))/𝐶_𝑝𝑛(true). 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.

(a)

(b)

(c)

(d)

Figure 4 

Comparison between the estimated 𝑘1, 𝑘2, 𝑘3, and 𝑘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 𝐶𝑝𝑛 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.

Figure 5 

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

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

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Comparison between the estimated 𝑘1, 𝑘2, 𝑘3, and 𝑘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 the initial peak is erroneous, and (e,f) when the tail part is erroneous.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Comparison between the estimated 𝑘1, 𝑘2, 𝑘3, and 𝑘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 initial peak is erroneous, and the bottom two when the tail part is erroneous.

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.

Table 1 

The percentage error in estimated 𝑘₁, 𝑘₂, 𝑘₃, and 𝑘₄ 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:(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.

Table 2 

The percentage error in estimated 𝑘₁, 𝑘₂, 𝑘₃, and 𝑘₄ 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.

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

Appendix

For a solution 𝑎,𝑏,𝑐,𝑑 that minimizes the cost function the partial derivatives of cost function with respect to the kinetic parameters is zero 𝜕0=𝑅̂̂𝑑,𝐶𝜕𝑎𝑎,𝑏,𝑐,_𝑝𝑎=𝑎,𝜕0=𝑅̂̂𝑑,𝐶𝜕𝑏𝑎,𝑏,𝑐,_𝑝𝑏=𝑏̂,𝜕0=𝑅̂̂𝑑,𝐶𝜕𝑐𝑎,𝑏,𝑐,_𝑝𝑐=𝑐̂,𝜕0=𝑅̂̂𝑑,𝐶𝜕𝑑𝑎,𝑏,𝑐,_𝑝𝑑=𝑑. (A.1) Let us define an implicit function which will be convenient for the application of chain rule 𝑎,𝑏̂,𝑐,̂𝑑=𝑔𝐶_𝑝=𝑔₁𝐶_𝑝,𝑔₂𝐶_𝑝,𝑔₃𝐶_𝑝,𝑔₄𝐶_𝑝^′ (A.2) that maps the 𝐶_𝑝 into an estimate [𝑎,𝑏̂,𝑐,̂𝑑]. We can rewrite (13) as 𝜕0=𝑅𝑔𝐶𝜕𝑎^𝑝,𝐶_𝑝,𝜕0=𝑅𝑔𝐶𝜕𝑏^𝑝,𝐶_𝑝,𝜕0=𝑅𝑔𝐶𝜕𝑐^𝑝,𝐶_𝑝,𝜕0=𝑅𝑔𝐶𝜕𝑑^𝑝,𝐶_𝑝. (A.3) By applying the chain rule to this equation, we can differentiate the above two equations with respect to 𝐶_𝑝 and obtain 𝜕0=²𝜕𝑎²𝑅𝑔𝐶_𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₁𝐶_𝑝,+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑏^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₂𝐶_𝑝,+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑐^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₃𝐶_𝑝,+𝜕²𝑅𝑔𝐶𝜕𝑎𝜕𝑑^𝑝,𝐶_𝑝𝜕𝜕𝐶_𝑝𝑔₄𝐶_𝑝,+𝜕²𝜕𝑎𝜕𝐶_𝑝𝑅𝑔𝐶_𝑝,𝐶_𝑝, (A.4) and similar equations for 𝑏, 𝑐, and 𝑑 can be obtained easily.

We can write these four equations in matrix form −𝜕²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅=𝜕²𝜕𝑎²𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅𝜕𝜕𝑎𝜕𝑑₂𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅𝜕𝜕𝑏𝜕𝑑₂𝑅𝜕𝜕𝑐𝜕𝑎²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝑅𝜕𝜕𝑑𝜕𝑎²𝑅𝜕𝜕𝑑𝜕𝑏²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅×𝜕𝑔₁𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₂𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₃𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₄𝐶_𝑝𝜕𝐶_𝑝𝑛. (A.5) Then a simple matrix inversion provides𝜕𝑔₁𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₂𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₃𝐶_𝑝𝜕𝐶_𝑝𝑛𝜕𝑔₄𝐶_𝑝𝜕𝐶_𝑝𝑛=𝜕²𝜕𝑎²𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅𝜕𝜕𝑎𝜕𝑑₂𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅𝜕𝜕𝑏𝜕𝑑₂𝑅𝜕𝜕𝑐𝜕𝑎²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝑅𝜕𝜕𝑑𝜕𝑎²𝑅𝜕𝜕𝑑𝜕𝑏²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅⁻¹×−𝜕²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅−𝜕²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅. (A.6) The elements of the matrix and the vector on the right hand side can be calculated based on the cost function 𝜕²𝑅𝜕𝑎²=2𝐶_𝑝^′𝜕𝐻𝜕𝑎^′𝜕𝐻𝐶𝜕𝑎^𝑝,𝜕²𝑅𝐶𝜕𝑎𝜕𝑏=2_𝑝^′𝜕𝐻𝜕𝑎^′𝜕𝐻𝐶𝜕𝑏^𝑝,𝜕²𝑅𝐶𝜕𝑎𝜕𝑐=−2_𝑝^′𝜕²𝐻𝜕𝑎𝜕𝑐^′𝐹+2𝐶_𝑝^′𝜕²𝐻𝜕𝑎𝜕𝑐^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑎^′𝜕𝐻𝐶𝜕𝑐^𝑝,𝜕²𝑅𝐶𝜕𝑎𝜕𝑑=2_𝑝^′𝜕𝐻𝜕𝑎^′𝜕𝐻𝐶𝜕𝑑^𝑝,𝜕²𝑅𝜕𝑏²=2𝐶_𝑝^′𝜕𝐻𝜕𝑏^′𝜕𝐻𝐶𝜕𝑏^𝑝,𝜕²𝑅𝐶𝜕𝑏𝜕𝑐=2_𝑝^′𝜕𝐻𝜕𝑏^′𝜕𝐻𝐶𝜕𝑐^𝑝,𝜕²𝑅𝐶𝜕𝑏𝜕𝑑=−2_𝑝^′𝜕²𝐻𝜕𝑏𝜕𝑑^′𝐹+2𝐶_𝑝^′𝜕²𝐻𝜕𝑏𝜕𝑑^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑏𝜕𝐻𝐶𝜕𝑑^𝑝,𝜕²𝑅𝜕𝑐²=−2𝐶_𝑝^′𝜕²𝐻𝜕𝑐^2′𝐹+2𝐶_𝑝^′𝜕²𝐻𝜕𝑐^2′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑐^′𝜕𝐻𝐶𝜕𝑐^𝑝,𝜕²𝑅𝐶𝜕𝑐𝜕𝑑=2_𝑝^′𝜕𝐻𝜕𝑐^′𝜕𝐻𝐶𝜕𝑑^𝑝,𝜕²𝑅𝜕𝑑²=−2𝐶_𝑝^′𝜕²𝐻𝜕𝑑^2′𝐹+2𝐶_𝑝^′𝜕²𝐻𝜕𝑑^2′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑑^′𝜕𝐻𝐶𝜕𝑑^𝑝,𝜕²𝑅𝜕𝑎𝜕𝐶_𝑝𝑛=−20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑎^′𝐹+20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑎^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑎^′𝐻0⋯01_𝑖=𝑛0⋯0^′,𝜕²𝑅𝜕𝑏𝜕𝐶_𝑝𝑛=−20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑏^′𝐹+20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑏^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑏^′𝐻0⋯01_𝑖=𝑛0⋯0^′,𝜕²𝑅𝜕𝑐𝜕𝐶_𝑝𝑛=−20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑐^′𝐹+20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑐^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑐^′𝐻0⋯01_𝑖=𝑛0⋯0^′,𝜕²𝑅𝜕𝑑𝜕𝐶_𝑝𝑛=−20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑑^′𝐹+20⋯01_𝑖=𝑛0⋯0𝜕𝐻𝜕𝑑^′𝐻𝐶_𝑝+2𝐶_𝑝^′𝜕𝐻𝜕𝑑^′𝐻0⋯01_𝑖=𝑛0⋯0^′, (A.7) where “′” denotes the transpose. The terms 𝜕𝐻/𝜕𝑎, 𝜕𝐻/𝜕𝑏,𝜕𝐻/𝜕𝑐, 𝜕𝐻/𝜕𝑑, (𝜕²/𝜕𝑎𝜕𝑐)𝐻, (𝜕²/𝜕𝑐²)𝐻, (𝜕²/𝜕𝑏𝜕𝑑)𝐻, and (𝜕²/𝜕𝑑²)𝐻 can be calculated using the compartment model 𝜕𝐻𝑒𝜕𝑎=T^−𝑐𝑡,𝜕𝐻𝑒𝜕𝑏=T^−𝑑𝑡,𝜕𝐻𝜕𝑐=T−𝑎𝑡𝑒^−𝑐𝑡,𝜕𝐻𝜕𝑑=T−𝑏𝑡𝑒^−𝑑𝑡,𝜕²𝜕𝑎𝜕𝑐𝐻=T−𝑡𝑒^−𝑐𝑡,𝜕²𝜕𝑏𝜕𝑑𝐻=T−𝑡𝑒^−𝑑𝑡,𝜕²𝜕𝑐²𝐻=T𝑎𝑡²𝑒^−𝑐𝑡,𝜕²𝜕𝑑²𝐻=T𝑏𝑡²𝑒^−𝑑𝑡, (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. 𝑠_1,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛,𝑠,𝑎,𝑏,𝑐,𝑑_2,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ2𝑠,𝑎+^1,12𝑠,𝑏+^1,22𝑠,𝑐+^1,32𝑠,𝑑+^1,42,𝑠_3,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ2𝑠,𝑎+^2,12𝑠,𝑏+^2,22𝑠,𝑐+^2,32𝑠,𝑑+^2,42,𝑠_4,𝑖=ℎ𝑢_𝑖𝐶_𝑝𝑛+ℎ,𝑎+𝑠_3,1,𝑏+𝑠_3,2,𝑐+𝑠_3,3,𝑑+𝑠_3,4,𝑦_𝑖𝐶_𝑝𝑛+ℎ=𝑦_𝑖𝐶_𝑝𝑛+𝑠_1,𝑖6+𝑠_2,𝑖3+𝑠_3,𝑖3+𝑠_4,𝑖6, (A.9) where 𝑘_𝑖=𝑦_𝑖(𝐶_𝑝𝑛), 𝑖=1,2,3,4 for 𝑎, 𝑏, 𝑐, and 𝑑, ℎ is the small step size we define according to the demand on accuracy and speed, 𝑢₁=1𝜙−𝜕²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅−𝜕𝜕𝑎𝜕𝑑₂𝜕𝑏𝜕𝐶_𝑝𝑛𝑅𝜕²𝜕𝑏²𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅−𝜕𝜕𝑏𝜕𝑑₂𝜕𝑐𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝑅𝜕²𝑅−𝜕𝜕𝑐𝜕𝑑₂𝜕𝑑𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑑𝜕𝑏²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅,𝑢₂=1𝜙𝜕²𝜕𝑎²𝜕𝑅−²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅𝜕𝜕𝑎𝜕𝑑₂𝜕𝜕𝑏𝜕𝑎𝑅−²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅𝜕𝜕𝑏𝜕𝑑₂𝜕𝜕𝑐𝜕𝑎𝑅−²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅𝜕²𝜕𝑐²𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝜕𝜕𝑑𝜕𝑎𝑅−²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅,𝑢₃=1𝜙𝜕²𝜕𝑎²𝑅𝜕²𝜕𝜕𝑎𝜕𝑏𝑅−²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑑₂𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝜕𝑅−²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑑₂𝑅𝜕𝜕𝑐𝜕𝑎²𝜕𝜕𝑐𝜕𝑏𝑅−²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝑅𝜕𝜕𝑑𝜕𝑎²𝜕𝜕𝑑𝜕𝑏𝑅−²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅𝜕²𝜕𝑑²𝑅,𝑢₄=1𝜙𝜕²𝜕𝑎²𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝜕𝜕𝑎𝜕𝑐𝑅−²𝜕𝑎𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝑅𝜕²𝜕𝜕𝑏𝜕𝑐𝑅−²𝜕𝑏𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑎²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝜕𝑅−²𝜕𝑐𝜕𝐶_𝑝𝑛𝑅𝜕²𝑅𝜕𝜕𝑑𝜕𝑎²𝑅𝜕𝜕𝑑𝜕𝑏²𝜕𝜕𝑑𝜕𝑐𝑅−²𝜕𝑑𝜕𝐶_𝑝𝑛𝑅, (A.10) and 𝜙 is a scalar given by the determinant 𝜕𝜙=²𝜕𝑎²𝑅𝜕²𝑅𝜕𝜕𝑎𝜕𝑏²𝑅𝜕𝜕𝑎𝜕𝑐²𝑅𝜕𝜕𝑎𝜕𝑑₂𝑅𝜕𝜕𝑏𝜕𝑎𝜕𝑏²²𝑅𝜕²𝑅𝜕𝜕𝑏𝜕𝑐²𝑅𝜕𝜕𝑏𝜕𝑑₂𝑅𝜕𝜕𝑐𝜕𝑎²𝑅𝜕𝜕𝑐𝜕𝑏𝜕𝑐²²𝑅𝜕²𝑅𝜕𝜕𝑐𝜕𝑑₂𝑅𝜕𝜕𝑑𝜕𝑎²𝑅𝜕𝜕𝑑𝜕𝑏²𝑅𝜕𝜕𝑑𝜕𝑐𝜕𝑑²²𝑅, (A.11)

Let us define 𝑎_𝑖=𝑦₁𝐶_𝑝𝑛,𝑏_𝑖=𝑦₂𝐶_𝑝𝑛,𝑐_𝑖=𝑦₃𝐶_𝑝𝑛,𝑑_𝑖=𝑦₄𝐶_𝑝𝑛,𝑎_𝑖+1=𝑦₁𝐶_𝑝𝑛,𝑏+ℎ_𝑖+1=𝑦₂𝐶_𝑝𝑛,𝑐+ℎ_𝑖+1=𝑦₃𝐶_𝑝𝑛,𝑑+ℎ_𝑖+1=𝑦₄𝐶_𝑝𝑛+ℎ, (A.12) where ℎ is as defined before. In multidimension, we perform the calculations 𝑁 times sequentially for every 𝐶_𝑝𝑛 for a single step 𝑎_𝑖+1=𝑎_𝑖+^𝑁_𝑛=1𝜕𝑎𝜕𝐶_𝑝𝑛Δ𝐶_𝑝𝑛,𝑏_𝑖+1=𝑏_𝑖+^𝑁_𝑛=1𝜕𝑏𝜕𝐶_𝑝𝑛Δ𝐶_𝑝𝑛,𝑐_𝑖+1=𝑐_𝑖+^𝑁_𝑛=1𝜕𝑐𝜕𝐶_𝑝𝑛Δ𝐶_𝑝𝑛,𝑑_𝑖+1=𝑑_𝑖+^𝑁_𝑛=1𝜕𝑑𝜕𝐶_𝑝𝑛Δ𝐶_𝑝𝑛. (A.13)