Cumulative Input Function Method for Linear Compartmental Models and Spectral Analysis in PET

Abstract

Compartmental modeling and spectral analysis are often used for tracer kinetic modeling in positron emission tomography (PET). The concentrations in kinetic equations are usually considered to be instantaneous, whereas PET data are inherently integrated over time, which leads to uncertainties in the results. A new formalism for kinetic analysis that uses cumulative tracer concentrations and avoids approximating the image-derived input function and PET measurements with midframe instantanous values was developed. We assessed the improvements of the new formalism over the midframe approximation methods for three commonly used radiopharmaceuticals: [¹¹C]raclopride, 2′-deoxy-2′-[¹⁸F]fluoro-d-glucose (FDG), and 3′-deoxy-3′-[¹⁸F]fluoro-thymidine (FLT). We found that improvements are case dependent and often not negligible. Improvements for determination of binding potential for [¹¹C]raclopride ranged from 5% to 25%. Improvements in estimation accuracy of FDG and FLT microparameters ranged up to 25%. On the other hand, estimation of macroparameter K _i = K₁k₃/(k₂ + k₃) for FDG or FLT did not show significant benefit with the new method; only modest improvement up to 2% was observed. Assessment of the benefits of using new method is far from being exhaustive, but possibly significant improvement was demonstrated. Therefore, we consider the proposed algorithm a necessary component of any kinetic analysis software.

Introduction

Compartmental models are frequently used to analyze dynamic positron emission tomography (PET) images for the kinetics of radiotracer tissue uptake. There exists several well-established linear compartmental models for the quantification of blood flow (Kety, 1951), neuroreceptor ligand binding (Lammertsma et al, 1996; Mintun et al, 1984), metabolic rate of glucose (Huang et al, 1980; Phelps et al, 1979; Sokoloff et al, 1977), and thymidine analog uptake and phosphorylation (Muzi et al, 2005a, b, 2006; Shields et al, 2002). Similarly, spectral analysis can also be used for the same purpose (Cunningham and Jones, 1993; Murase, 2003). Although it differs from compartmental modeling conceptually, it does not differ much technically, because in both cases, the solution for the tissue time activity curve (tTAC) consists of a linear combination of exponential functions and possibly a constant term, convolved with the input function (IF) (Gunn et al, 2001; Schmidt, 1999).

Unlike conventional physiological sampling, PET scanning does not allow measurement of instantaneous concentrations, because it measures them indirectly, through imaging tracer activity over a predefined time period. That issue was fully solved for the graphical methods (Buchert et al, 2003), but no general solution for compartmental models and spectral analysis has been proposed. A popular kinetic analysis software package, SAAM II (Barrett et al, 1998; Berneis et al, 2004) simply ignores that issue and assumes all concentration measurements as instantaneous. That approximation is often called midframe approximation; that is, the TAC in the middle of the time frame is equalized to the TAC averaged over that time frame. Other popular kinetic analysis software packages, including COMKAT (Muzic and Cornelius, 2001) and KMZ (Burger and Buck, 1997), now called PMOD (PMOD Technologies Ltd., Zürich, Switzerland), explicitly integrate instantaneous model output over the corresponding time period, which is appropriate for the blood sampling IF but does not represent the PET image-derived IF correctly.

We have developed a cumulative concentration IF approach that is especially suited for the kinetic models with cumulative IF measurements like PET, but also useful for instantaneous measurements. We have investigated the potential improvements of the parameter estimation accuracy compared with the midframe approximations for several well-established PET compartmental models.

Theoretical Background

The general linear compartmental system with n compartments is described by the state equation (equation 1).

$$\begin{array}{rl}& \dot{\mathbf{c}}(t)=\mathbf{A} \mathbf{i}+\mathbf{K} \mathbf{c}(t)\\ & \mathbf{y}(t)=\mathbf{B} \mathbf{c}(t)+\mathbf{D} \mathbf{i}(t)\\ & \mathbf{c}(0)={\mathbf{c}}_{\mathbf{0}}\end{array}$$

(1)

In equation 1, c(t) is an n dimensional column vector of state variables (i.e., concentrations in the compartments), i(t) is a p-dimensional column vector of IFs and y(t) is a q-dimensional column vector of observed outputs (Jacquez, 1996). A is an n × p matrix of influx rate parameters from each input to each compartment. The state transition matrix K is an n × n matrix that consists of rate parameters k _ij , which describe transfer from the compartment j to the compartment i when j ≠ i. The k _ii stands for the sum of all the transfers from the compartment i to other compartments (k_ji in equation 2), for the transfer from the compartment i to the environment (k _i in equation 2), and for the decrease of specific activity of the radiotracer because of the nuclear decay (λ in equation 2), if the data are not decay corrected.

$$k_{ii}=-\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{k}_{ji}-k_i-\lambda$$

(2)

B is a q × n observation matrix and D is a q × p feedthrough matrix that describes the direct contribution of each input to the observed output. The elements of A, K, B, and D in equation 1 are assumed to be constant with respect to time during the imaging period; their values may change from one experimental situation to another. The solution of the state equation can be described as

$$\begin{array}{rl}& \mathbf{c}(t)=\underset{0}{\overset{\mathrm{\infty }}{\int }}{\mathrm{e}}^{\mathbf{K}(\tau )}\mathbf{A}\mathbf{i}(t-\tau )d\tau \\ & \mathbf{y}(t)=\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathbf{B}{\mathrm{e}}^{\mathbf{K}(\tau )}\mathbf{A}\mathbf{i}(t-\tau )d\tau +\mathbf{D} \mathbf{i}(t)\\ & {\mathrm{e}}^{\mathbf{K}\tau }=\mathbf{J}+\tau \mathbf{K}+({\tau }^2/2!){\mathbf{K}}^2+({\tau }^3/3!){\mathbf{K}}^3+\cdots \end{array}$$

(3)

The term e^Kτ in equation 3 is the matrix exponential and J is the n × n identity matrix. As we are interested in time integral of measurements over some predefined time interval [t_j–1, t _j ], we integrate the second equation on equation 3.

$$\begin{array}{rl}\mathbf{Y}(t)& =\underset{-\mathrm{\infty }}{\overset{t}{\int }}\mathbf{y}(\tau )\mathrm{d}\tau \\ & =\underset{-\mathrm{\infty }}{\overset{t}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathbf{B}{\mathrm{e}}^{\mathbf{K}(\xi )}\mathbf{A}\mathbf{i}(\tau -\xi )\mathrm{d}\xi \mathrm{d}\tau +\underset{-\mathrm{\infty }}{\overset{\mathrm{t}}{\int }}\mathbf{D} \mathbf{i}(\tau )\mathrm{d}\tau \end{array}$$

(4)

In equation 5, we exchange the order of integration. By rewriting equation 4 and imposing cumulative IF vector I(t) we get:

$$\begin{array}{rl}& \underset{-\mathrm{\infty }}{\overset{t}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathbf{B}{\mathrm{e}}^{\mathbf{K}(\xi )}\mathbf{A}\mathbf{i}(\tau -\xi )\mathrm{d}\xi \mathrm{d}\tau =\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathbf{B}{\mathrm{e}}^{\mathbf{K}(\xi )}\mathbf{A}\underset{-\mathrm{\infty }}{\overset{t-\xi }{\int }}\mathbf{i}(u) \mathrm{d}u\mathrm{d}\xi \\ & \mathbf{I}(t)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}\mathbf{i}(\tau )\mathrm{d}\tau \\ & \mathbf{Y}(t)=\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathbf{B}{\mathrm{e}}^{\mathbf{K}(\xi )}\mathbf{A}\mathbf{I}(t-\xi )\mathrm{d}\xi +\mathbf{D}\mathbf{I}(t)\end{array}$$

(5)

Comparing equation 5 to equation 3, we notice that the relation for cumulative vector of observed outputs Y(t) is almost the same as the relation for instantaneous vector of observed outputs y(t), except that the IF in the first case is cumulative, whereas it is instantaneous in the second case. Time average of instantaneous vector of observed output y(t) over some predefined time interval [t_j–1, t _j ] is easily calculated using cumulative vector of observed outputs Y(t)

$$\mathbf{z}(t_{j-1},t_j)=\frac{1}{t_j-t_{j-1}}\underset{t_{j-1}}{\overset{t_j}{\int }}\mathbf{y}(\tau )\mathrm{d}\tau =\frac{\mathbf{Y}(t_{j-1})-\mathbf{Y}(t_j)}{t_j-t_{j-1}}$$

(6)

Positron Emission Tomography Compartmental Modeling

In the PET compartmental modeling case, there is one observed output and one or more inputs. The observed output z(t_j–1, t _j ) is always time averaged over time frame [t_j–1, t _j ], but a kinetic model can easily take this into account by explicit integration of the instantaneous model output. Input functions are instantaneous in the case of blood sampling and time averaged over time frames when they are image derived. Exact conversion from instantaneous to time-averaged IF is possible, whereas the opposite is impossible. Depending on the available experimental data, we can apply three types of midframe approximations:

A. correct IF, approximated tTAC,

B. approximated IF, correct tTAC,

C. approximated IF and approximated tTAC.

Improvement of parameter estimation accuracy when using the new method instead of the A, B, or C approximations could be quantified by finding most probable parameter estimations for the approximate methods and calculating the difference between these estimations and correct parameters. Most probable parameter estimations for the approximate methods are found by minimizing the differences between correct model output representation z(t_j–1, t _j , θ), with correct parameters θ and approximate model output representation Z(t_j–1, t _j , $\hat{\theta }$ ), with most probable parameters $\hat{\theta }$ . Nonlinear problem with objective function Φ in equation 7 is solved by one of the standard methods for minimum search of multivariable function.

$$\mathrm{\Phi }(\hat{\theta })=\sum _i{w}_i(\hat{z}(t_{i-1},t_i,\hat{\theta })-z(t_{i-1},t_i,\theta ){)}^2$$

(7)

Weights w _i in equation 7 must be the same as we would use for parameter estimation with weighted least squares method.

Materials and methods

The improvement in parameter estimation with the new method, compared with all three possible approximations A, B, and C, was quantified for three common radiopharmaceuticals; [¹¹C]raclopride, 2′-deoxy-2′-[¹⁸F]fluoro-d-glucose (FDG), and the 3′-deoxy-3′-[¹⁸F]fluoro-thymidine (FLT). Input function and scanning protocol was obtained from the reference study, listed further in text for each case. In our analysis, the IF was represented with splines of fourth order and considered as errorless representation. Such representation enabled sampling of the instantaneous IF, its cumulative value and its average over time frames (midframe approximation). The expression for the model output, which is calculated as a convolution between IF and model kernel, was derived analytically and evaluated for the new method and all three approximations. We set up a nonlinear system (equation 7) for one typical parameter value set for each radiopharmaceutical (listed in Tables 1 to 3) and solved it with the weights, calculated as a square root of the ratio between time frame duration and model output that has not been decay corrected. All the calculations were performed in Matlab (The MathWorks, Natick, MA, USA).

Table 1

Accuracy improvements for the [¹¹C]raclopride

Parameter	K1 (mL/min/g)	k2 (1/min)	k3 (1/min)	BPND = k3/k4	V b	V d
Correct value	0.1	0.11	0.12	1.1	0.05	1.91
A
Estimation	0.099	0.107	0.111	1.05	0.049	1.91
Δabsolute	−0.0006	−0.0033	−0.0086	−0.049	−0.0013	0.0008
Δrelative	−0.59%	−3.0%	−7.2%	−4.5%	−2.6%	0.04%
B
Estimation	0.094	0.092	0.096	0.86	0.052	1.90
Δabsolute	−0.0056	−0.018	−0.024	−0.24	0.0022	−0.0056
Δrelative	−5.6%	−16%	−20%	−22%	4.4%	−0.29%
C
Estimation	0.094	0.092	0.093	0.86	0.051	1.90
Δabsolute	−0.0056	−0.018	−0.027	−0.24	0.0014	−0.0051
Δrelative	−5.6%	−16%	−22%	−22%	2.8%	−0.27%

Table shows the correct parameter values, the parameters we would estimate with all three possible approximations, according to our methodology and the absolute and relative improvements with new method. BP_ND is the most relevant clinical parameter and our method demonstrates the improvement in estimation accuracy in the range of 5% to 25%. Alternative calculation of this parameter uses the volume of distribution in striatum and cerebellum. Therefore, we list also the improvement in estimation accuracy for volume of distribution (V_d), which is well below 1%.

Table 2

Accuracy improvements for the FDG

Parameter	K1 (mL/min/g)	k2 (1/min)	k3 (1/min)	k4 (1/min)	V b	K i (mL/min/g)
Correct value	0.07	0.09	0.05	0.005	0.02	0.025
A
Estimation	0.070	0.090	0.051	0.0054	0.020	0.025
Δabsolute	0.00001	0.0004	0.0008	0.0004	−0.0004	0.0002
Δrelative	0.01%	0.5%	1.5%	7.7%	−2.2%	0.7%
B
Estimation	0.072	0.100	0.053	0.0056	0.015	0.025
Δabsolute	0.002	0.007	0.003	0.0006	−0.005	0.0004
Δrelative	2.9%	7.3%	5.3%	12%	−27%	1.7%
C
Estimation	0.070	0.090	0.050	0.0051	0.019	0.025
Δabsolute	−0.0002	−0.0001	0.0001	0.0001	−0.001	−0.0000
Δrelative	−0.22%	−0.10%	0.24%	2.6%	−5.7%	−0.006%

While the proposed method does not substantially change macroparameter K _i estimation, microparameters estimation such as V_b may show more substantial change.

Table 3

Accuracy improvements for the FLT

Parameter	K1 (mL/min/g)	k2 (1/min)	k3 (1/min)	k4 (1/min)	V b	K i (mL/min/g)
Correct value	0.1	0.1	0.1	0.02	0.05	0.05
A
Estimation	0.100	0.098	0.098	0.020	0.049	0.050
Δabsolute	−0.0005	−0.0023	−0.0022	−0.0001	−0.0013	−0.0002
Δrelative	−0.47%	−2.3%	−2.2%	−0.39%	−2.6%	−0.42%
B
Estimation	0.100	0.101	0.099	0.020	0.050	0.050
Δabsolute	−0.00005	0.0007	−0.0005	−0.0002	0.00005	−0.0003
Δrelative	−0.050%	0.69%	−0.53%	−1.2%	0.093%	−0.66%
C
Estimation	0.100	0.101	0.100	0.020	0.049	0.050
Δabsolute	0.0001	0.001	−0.0004	−0.0003	−0.0005	−0.0003
Δrelative	0.07%	1.1%	−0.44%	−1.5%	−1.0%	−0.67%

Like the FDG case, the proposed method does not substantially change macroparameter K _i estimation. However, also the microparameter estimation changes less than they change in the FDG case, although still much more than the macroparameter K _i estimation.

For [¹¹C]raclopride, the improvement of the new method was estimated for the two-compartment four-parameter model with the plasma IF and the following kinetic parameters K₁, k₂, k₃, BP_ND = k₃/k₄, and V_b. This model is used for the estimation of the binding potential, a measure of f_NDB_avail/K_d, of [¹¹C]raclopride in human striatum. Alternative calculation of binding potential uses the volume of distribution V_d = k₁/k₂(1 + k₃/k₄) in striatum and cerebellum. Therefore, we calculate also the improvement in estimation accuracy for the volume of distribution. We obtained IF and striatal TAC (shown in Figure 1), as well as scanning protocol from Lammertsma et al (1996).

Figure 1

Input function and tissue time activity curve (tTAC) for the [¹¹C]raclopride simulation.

For FDG and FLT, the improvement of the new method was estimated for the two-compartment four-parameter kinetic model with the plasma IF and the following microparameters K₁, k₂, k₃, k₄, and V_b. Improvement in estimation of macroparameter K _i = K₁k₃/(k₂ + k₃), which defines metabolic rate of glucose in case of FDG PET and thymidine blood-tissue transfer constant in case of FLT PET, was also quantified. We obtained IF and scanning protocol for FDG from Schroeder et al (2007) and for FLT from Muzi et al (2005a).

Results

Figure 2 shows the differences between the correct tTAC and approximate tTACs that are calculated with correct parameters or with the most probable approximate method parameter estimations, according to equation 7 for all three approximations A, B, and C. Differences between the correct tTAC, calculated with correct parameters, and approximate tTACs, calculated with the most probable approximate method parameter estimation, are considerable only shortly after injection, when IF and tTAC changes rapidly, and only for approximation A. These differences are lower for the approximation B and almost negligible for the approximation C. They are also negligible for all three approximations at a later time after injection. Differences between the correct tTAC, calculated with correct parameters, and approximate tTACs, calculated with the most probable approximate method parameter estimation, are also considerable for approximation A shortly after injection. They are lower for the approximation B and almost negligible for the approximation C. At later time after injection, these differences are the highest for the approximation C and lowest for the approximation A. They are considerably higher than the differences between the correct tTAC, calculated with correct parameters, and approximate tTACs, calculated with the most probable approximate method parameter estimation. This observation implies that minimizing objective function in equation 7 was successful and alignment between correct model output and approximate model output was improved.

Figure 2

The assessment of the new method improvements—data for the [¹¹C]raclopride. Graphs show the differences between correct tissue time activity curve (tTAC) and the tTACs that are calculated with one of the approximation methods, either with the correct parameters (dashed line) or with the parameters that maximizes alignment with the correct tTAC (solid line). (A) Approximation A. (B) Approximation B. (C) Approximation C.

Table 1 shows the results for [¹¹C]raclopride simulation. Absolute and relative differences are recognized as an improvement of new method, compared with approximations. Improvements in microparameter estimation accuracy are moderate for approximation A, ranging from 1% for K₁ to 8% for k₃, whereas improvement in volume of distribution estimation accuracy is only 0.04%. Estimation accuracy for the most relevant clinical parameter BP_ND = k₃/k₄ is improved for 5%. Considerable higher improvements are for the methods B and C. Improvements are well above 1% for all parameters, whereas BP_ND= k₃/k₄ estimation accuracy is improved for 20% to 25%. Alternative calculation of this parameter is by using the formula BP_ND = (V_t–V_ref)/V_ref, where V_t is the volume of distribution in striatum and V_ref is the volume of distribution in cerebellum. Improvement in estimation accuracy for volume of distribution is well below 1% also for methods B and C.

Tables 2 and 3 show the results for FDG and FLT simulation. Unlike the [¹¹C]raclopride case, where the approximation A was more accurate, for FDG and FLT the accuracy of all three approximations is similar. One exception might be approximation B for the FDG case, which is noticeably less accurate than others. Estimation of macroparameter K _i is fairly accurate with all three approximations. New method improves estimation accuracy above 1% only in the approximation B for FDG case. Estimation of microparameters is slightly less accurate and new method improves parameter estimation accuracy often above 1%. Improvements in V_b estimation reaches 27% in the approximation B for FDG case.

Discussion

A brief check of the results in tables reveals that significance of the new method improvement is case dependent. The [¹¹C]raclopride case shows that the new method can lead to significant parameter estimation accuracy improvement, compared with the midframe approximations. The [¹¹C]raclopride is used for the estimation of the binding potential, a measure of f_NDB_avail/K_d, in human striatum. Therefore, the most interesting kinetic parameter is BP_ND = k₃/k₄. It is evident from Table 1 that the improvement in BP_ND estimation accuracy when using new method was in the range of 5% to 25%. Alignment plots in Figure 2 show that the most probable approximate method parameter estimation improves the alignment between approximate and correct tTAC, although the alignment is still not perfect. Therefore, image noise and other sources of errors could also have the significant effect on the estimation of binding potential in human striatum. Another scanning protocol, IF shape, and parameter values set would probably lead to different values, but this example alone indicates that potential improvement in BP_ND estimation accuracy when using our new method can be very high.

Binding potential can be calculated also by using the formula BP_ND = (V_t—V_ref)/V_ref. It is interesting that the improvement in volume of distribution estimation accuracy is almost negligible (well below 1%), which indicates that improvement in alternative binding potential estimation accuracy is small.

In FDG and FLT kinetic analysis, we are usually interested in parameter K _i , which is a macroparameter and generally much more stable than microparameters. We found that the new method improves its estimation accuracy only slightly—below 2% in our case. Parameter estimation accuracy improvements for other FDG or FLT cases with different scanning protocols, IF shapes, and parameter values set might differ from our case. However, because FDG or FLT uptake dynamics are quite linear after a certain time, the midframe approximation is expected not to have a major effect. On the other hand, as shown in the case [¹¹C]raclopride, the further from linear the uptake dynamics is, the more significant effects one can expect. From the results in Tables 2 and 3, it is also evident that the possible improvement in microparameter estimation accuracy is > 25%, which is much higher than the K _i estimation accuracy improvements. It is interesting that we got significantly different results for FDG and FLT case, although compartmental model is the same. Differences in results are attributed to different IF, scanning protocol, and simulated parameters.

We should be aware that these improvements in parameter estimation accuracy may be affected with some other factors, including model inadequacy and noise in PET measurements. These two factors impose some intrinsic parameter estimation inaccuracy, which interferes with the inaccuracy, arising from the midframe model approximations. As these sources of errors in parameter estimation method are uncorrelated, these errors add in quadrature, which decreases the effect of one source of error to the whole error.

Conclusion

One of three possible types of midframe approximation is typically applied in kinetic analysis, without being justified with error analysis or considered as a possible source of error. In this article, we proposed an algorithm that avoids the midframe approximation for cases when data are partially or fully in cumulative form. We also developed a method that quantifies the improvement in parameter estimation when avoiding midframe approximations, and applied it on three common radiopharmaceuticals; [¹¹C]raclopride, FDG, and FLT.

We found that the improvement in parameter estimation when avoiding midframe approximations is strongly case dependent and could be significant. If the midframe approximation is applied, it has to be justified. However, our newly developed algorithm is general, and it can be applied to various linear compartmental models with constant coefficients, with plasma IF or reference region IF, as well as to spectral analysis. The implementation of proposed algorithm does not increase computational complexity but rather decreases it. Therefore, we consider the proposed algorithm as a necessary improvement in implementation of any kinetic analysis software that relies on midframe approximations.