Multilinear models in PET data analysis

Non-linear fitting algorithms are used to estimate parameters of compartmental models (CM) from regional PET data. For pixel-by-pixel computation of parametric maps these methods tend to be too slow, and therefore much faster linear analysis methods are preferably used, including multiple-time graphical analysis (Logan and Patlak plots) and multilinear models.

The compartmental models for PET data are described in terms of a set of linear, first-order, constant-coefficient, ordinary differential equations (ODEs). In a reversible one-tissue compartment model (1TCM) one ODE is sufficient to describe the exchange between the tissue and arterial input compartments:

$\frac{dC_{1}(t)}{dt} = K_{1} C_{0}(t) - k_{2} C_{1}(t)$

Assuming that initial concentration in tissue compartment is zero (C₁(0)=0), Eq 1 can be integrated to provide the tissue concentration at time T:

$C_{1}(T) = K_{1} \int_{0}^{T} C_{0}(t)dt - k_{2} \int_{0}^{T} C_{1}(t)dt$

The radioactivity concentration in tissue, measured using PET, is contaminated by spillover from adjacent blood vessels and vascular volume inside the region-of-interest (or image voxel):

$C_{PET}(t) = V_{B} C_{B}(t) + C_{1}(t)$

Substituting C₁ with Eq 2 gives

$C_{PET}(t) = V_{B} C_{B}(t) + K_{1} \int_{0}^{T} C_{0}(t)dt - k_{2} \int_{0}^{T} C_{1}(t)dt$

Integration and rearrangement of Eq 3 gives

$\int_{0}^{T} C_{1}(t)dt = \int_{0}^{T} C_{PET}(t)dt - V_{B} \int_{0}^{T} C_{B}(t)dt$

which is substituted in Eq 4, providing equation

$C_{PET}(t) = V_{B} C_{B}(t) + K_{1} \int_{0}^{T} C_{0}(t)dt - k_{2} \int_{0}^{T} C_{PET}(t)dt + k_{2} V_{B} \int_{0}^{T} C_{B}(t)dt$

If the radioactivity concentration in blood can be represented by the model input function, as is the case with [¹⁵O]H₂O, then C_B=C₀, and the Eq 6 simplifies to Eq 7. This assumption can be made also with some other radiotracers, including [¹⁸F]FDG.

$C_{PET}(t) = V_{0} C_{0}(t) + \left( K_{1} + V_{0} k_{2} \right) \int_{0}^{T} C_{0}(t)dt - k_{2} \int_{0}^{T} C_{PET}(t)dt$

The three coefficients in Eq 7, p₁=V₀, p₂=K₁+V₀k₂, and p₃=k₂, can be estimated by multilinear regression

$y = p_1 x_1 + p_2 x_2 + p_3 x_3 \nonumber$

, where p_i denotes the coefficient of independent variable x_i. The number of coefficients and model parameters must be the same for the linearization method to be useful.

Since the p_i coefficients are either ≥0 or ≤0, the non-negative least squares (NNLS) method is well-suited for the estimation. As its input, NNLS requires matrix A and vector b, filled with y and x values, respectively. Values from time frames 1 - n are placed on their own rows:

$\begin{align*} \overset{A}{ \begin{bmatrix} C_{0}(t_1) & \int_{0}^{t_1} C_{0}(t)dt & -\int_{0}^{t_1} C_{PET}(t)dt \\ C_{0}(t_2) & \int_{0}^{t_2} C_{0}(t)dt & -\int_{0}^{t_2} C_{PET}(t)dt \\ \vdots & \vdots & \vdots \\ C_{0}(t_n) & \int_{0}^{t_n} C_{0}(t)dt & -\int_{0}^{t_n} C_{PET}(t)dt \end{bmatrix} } \overset{b}{ \begin{bmatrix} C_{PET}(t_1) \\ C_{PET}(t_2) \\ \vdots \\ C_{PET}(t_n) \end{bmatrix} } \end{align*}$

NNLS algorithm returns p_i coefficient vector, from which the compartment model parameters can be solved, in this case simply as:

$\begin{align*} V_0 &= p_1 \\ K_1 &= p_2 - p_1 p_3 \\ k_2 &= p_3 \end{align*}$

The multilinear solution method has been extended to two-tissue compartment model (2TCM) by Blomqvist (1984) and Evans (1987), and further for [¹⁸F]FDOPA model with rate constants K₁-k₅ (Gjedde, 1991). Equations for three-tissue compartmental models with six rate constants are given in TPCMOD0023 and TPCMOD0024. Equations for the irreversible and reversible multilinear 2TCMs are

$C_{PET}(t) = V_{0} C_{0}(t) + \left( K_1 + V_0 ( k_2 + k_3 ) \right) \int_{0}^{T} C_{0}(t)dt - \left( k_2 + k_3 \right) \int_{0}^{T} C_{PET}(t)dt + K_1 k_3 \int_{0}^{T} \int_{0}^{T} C_{0}(t)dt dt$

and

$C_{PET}(t) = V_0 C_{0}(t) + \left( K_1 + V_0 ( k_2 + k_3 + k_4 ) \right) \int_{0}^{T} C_{0}(t)dt - \left( k_2 + k_3 + k_4 \right) \int_{0}^{T} C_{PET}(t)dt \\ \quad \quad \quad \quad + \left( K_1 ( k_3 + k_4 ) + V_0 k_2 k_4 \right) \int_{0}^{T} \int_{0}^{T} C_{0}(t)dt dt - k_2 k_4 \int_{0}^{T} \int_{0}^{T} C_{PET}(t)dt dt$

The net influx rate K_i can be calculated from the parameters of the irreversible 2TCM model, and from the coefficient vector of the multilinear solution in Eq 8 as:

$K_i = \frac{K_1 k_3}{k_2 + k_3} = \frac{p_4}{p_3} \nonumber$

The distribution volume V_T can be calculated from the parameters of the reversible 2TCM model, and from the coefficient vector of the multilinear solution in Eq 9 as:

$V_T = \frac{K_1}{k_2} \left( 1 + \frac{k_3}{k_4} \right) = \frac{p_4}{p_3} - p_1 \nonumber$

Thus, the calculation of macroparameters K_i or V_T from the p_i coefficients or rate constants requires division, which leads to very high variation and too noisy parametric maps. Rearranging of equations can provide the macroparameters directly with much less noise although with somewhat increased biases (Zhou et al., 2003; Kim et al., 2008). For the irreversible 2TCM (Eq 8), the solution of rearranged equation

$\int_{0}^{T} C_{PET}(t)dt = \frac{V_0}{k_2 + k_3} C_{0}(t) - \frac{1}{k_2 + k_3} C_{PET}(t) + \left( \frac{K_1}{k_2 + k_3} + V_0 \right) \int_{0}^{T} C_{0}(t)dt \\ \quad \quad \quad \quad \quad \quad \quad + \frac{K_1 k_3}{k_2 + k_3} \int_{0}^{T} \int_{0}^{T} C_{0}(t)dt dt$

provides the K_i estimate directly as p₄. For the reversible 2TCM (Eq 9), the solution of rearranged equation

$\int_{0}^{T} \int_{0}^{T} C_{PET}(t)dt dt = \frac{V_0}{k_2 k_4} C_{0}(t) - \frac{1}{k_2 k_4} C_{PET}(t) + \left( \frac{K_1 + V_0 ( k_2 + k_3 + k_4 )}{k_2 k_4} \right) \int_{0}^{T} C_{0}(t)dt \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad - \left( \frac{k_2 + k_3 + k_4}{k_2 k_4} \right) \int_{0}^{T} C_{PET}(t)dt + \left( \frac{K_1}{k_2}(1 + \frac{k_3}{k_4}) \right) \int_{0}^{T} \int_{0}^{T} C_{0}(t)dt dt$

provides the V_T estimate without division as p₅ - p₁.

Reference region input compartmental models

Simplified reference tissue model (SRTM) can be linearized to compute BP_ND maps (Zhou et al., 2003):

$C_{T}(t) = R_1 C_{R}(t) + k_2 \int_{0}^{T} C_{R}(t)dt - \frac{k_2}{{BP}_{ND} + 1} \int_{0}^{T} C_{T}(t)dt$

, where R₁ is the ratio of K₁s in the region of interest and in the reference region, and by assumption of reference region models, also the ratio of k₂s:

$R_1 = \frac{K_1}{K_1^\prime} = \frac{k_2}{k_2^\prime}$

To reduce noise in BP_ND, the linearized model can be restructured to compute DVR (BP_ND+1) without division:

$\int_{0}^{T} C_{T}(t)dt = \left( {BP}_{ND} + 1 \right) \int_{0}^{T} C_{R}(t)dt + \frac{R_1}{k_2} \left( {BP}_{ND} + 1 \right) C_{R}(t) - \frac{1}{k_2} \left( {BP}_{ND} + 1 \right) C_{T}(t)$

Since the k₂' should not vary across brain pixels, in SRTM2 the SRTM model is constrained by fixing k₂' to its global median (Wu & Carson, 2002; Ichise et al., 2003). In multilinear version (MRTM2), the k₂' (=k₂/R₁) could first be estimated by minimizing model

$C_{R}(t) = \frac{1}{R_1} C_{T}(t) + \frac{k_2^\prime}{{BP}_{ND} + 1} \int_{0}^{T} C_{T}(t)dt - k_2^\prime \int_{0}^{T} C_{R}(t)dt$

and then, using the global median of k₂', BP_ND map can be computed by minimizing

$\int_{0}^{T} C_{T}(t)dt = \left( {BP}_{ND} + 1 \right) \left[ \frac{C_{R}(t)dt}{k_2^\prime} \right] - \left( \frac{{BP}_{ND} + 1}{R_1} \right) \frac{C_{T}(t)}{k_2^\prime}$

The transport-limited reference tissue model (TRTM), where the reference region is such that irreversible binding or metabolism is very rapid (k₃'>>k₂'), can be linearized to produce k₃ maps (TPCMOD0002):

$C_{T}(t) = R_1 C_{R}(t) + R_1 k_3 \int_{0}^{T} C_{R}(t)dt - \left( k_2 + k_3 \right) \int_{0}^{T} C_{T}(t)dt$

or to estimate k₃ directly without division:

$C_{R}(t) = \frac{1}{R_1} C_{T}(t) + \frac{k_2 + k_3}{R_1} \int_{0}^{T} C_{T}(t)dt - k_3 \int_{0}^{T} C_{R}(t)dt$

Applications

One-tissue compartment model (1TCM) is sufficient to describe the kinetics of many radiotracers, including radiowater. Linearization of radiowater model was described by van den Hoff et al (1993), and it was applied to skeletal muscle [¹⁵O]H₂O PET data by Burchert et al (1997), and to myocardial muscle by Lee et al (2005). Linearized 1TCM has been used to assess BBB permeability from [⁶⁸Ga]EDTA PET data (Zhou et al., 2001) and histamine receptor occupation with [¹¹C]dozepin PET (Zhou et al., 2002).

Multilinear form of irreversible 2TCM (k₄=0) and reversible 2TCM (k₄>0) has been used to compute parametric images from PET studies where 1TCM is not suitable, for example [¹¹C]glucose, [¹¹C]deoxyglucose, and [¹¹C]methionine (Blomqvist, 1984), [¹⁸F]FDG (Cai et al., 2002; Huang et al., 2007; Pouzot et al., 2013), [¹⁸F]fluoride (Kim et al., 2007; Sanchez-Crespo et al., 2017), dopamine 2 receptor radioligand [¹¹C]FLB457 (Hagelberg et al., 2004), and μ opioid receptor radioligand [¹¹C]MeNTI (Kim et al., 2008) PET data.

SRTM can be linearized to compute BP_ND maps (Zhou et al., 2003; TPCMOD0002). The method has been used in analyses of brain receptor studies, for example [¹¹C]DASB (Kim et al., 2006; Endres et al., 2011).

Time delay between arterial blood and tissue curves can be rapidly estimated using linearized model (van den Hoff et al., 1993).

Bias

Noise in the data leads to bias in parameter values estimated using linear least squares methods, (Feng et al., 1993 and 1996). The results of multilinear model fitting can be used as initial values for bias-free but slower non-linear optimization methods (Feng et al., 1993) or refined with general linear least squares (GLLS) method) (Feng & Ho, 1993; Feng et al., 1995 and 1996; Negoita & Renaut, 2005; Wen et al., 2007).

Multilinear equations with time frames

The radioactivity concentration curves (time-activity curves, TACs), measured with PET, consist of consecutive time frames; we do not have measurements of radioactivity concentrations at specific time points (as assumed in the equations above), but instead we have the average concentrations during each time frame, and the AUCs at the end of each time frame. At the time zero, AUC=0. Assuming that the frames are contiguous (the start time of frame i equals the end time of frame i-1), at the end of frame i the AUC is

${AUC}_{i}^{e} = {AUC}_{i-1}^{e} + {\Delta t}_i \times C_i \nonumber$

, where Δt is the length of time frame (Δt = t^e_i - t^e_i-1) and C_i is the average concentration during the frame i.

Regarding the one-tissue compartment model in Eq 2, the concentrations in the tissue compartment C₁ at the start and end of frame i are:

$C_{1,i}^{s} = C_{1,i-1}^{e} = K_{1} {AUC}_{0,i-1}^{e} - k_{2} {AUC}_{1,i-1}^{e} \nonumber$ $C_{1,i}^{e} = K_{1} {AUC}_{0,i}^{e} - k_{2} {AUC}_{1,i}^{e} \nonumber$

Substitution of Eq 5 into these gives

$C_{1,i}^{s} = ( K_{1} + V_0 k_2 ) {AUC}_{0,i-1}^{e} - k_{2} {AUC}_{PET,i-1}^{e} \nonumber$ $C_{1,i}^{e} = ( K_{1} + V_0 k_2 ) {AUC}_{0,i}^{e} - k_{2} {AUC}_{PET,i}^{e} \nonumber$

The average concentration during frame i is

$\begin{align*} C_{PET,i} &= V_{0} C_{0,i} + \frac{C_{1,i}^{s} + C_{1,i}^{e}}{2} \\ &= V_{0} C_{0,i} + ( K_{1} + V_0 k_2 ) \frac{ ({AUC}_{0,i-1}^{e} + {AUC}_{0,i}^{e}) }{2} \ - k_{2} \frac{ ({AUC}_{PET,i-1}^{e} + {AUC}_{PET,i}^{e}) }{2} \end{align*}$

The average AUC during frame i is

${AUC}_{i} = \frac{ {AUC}_{i-1}^{e} + {AUC}_{i}^{e} }{2} \nonumber$

and thus the PET time frame based multilinear equation for one-tissue compartment can be written as

$C_{PET,i} = V_{0} C_{0,i} + ( K_{1} + V_0 k_2 ) {AUC}_{0,i} - k_{2} {AUC}_{PET,i}$

, where all concentrations and AUCs represent the averages during the frame i. This equation has the same form as Eq 7 for the non-framed data. Equations for the other compartmental models can be similarly applied to PET data with time frames. The AUC average during PET frame i can be calculated as

${AUC}_{i} = {AUC}_{i-1}^{e} + \frac{\Delta t}{2} C_{i} \nonumber$

Literature

Blomqvist G. On the construction of functional maps in positron emission tomography. J Cereb Blood Flow Metab. 1984; 4:629-632. doi: 10.1038/jcbfm.1984.89.

Blomqvist G, Pauli S, Farde L, Eriksson L, Person A, Halldin C. Dynamic models for reversible ligand binding. In: Positron Emission Tomography in Clinical Research and Clinical Diagnosis: Tracer Modelling and Radioreceptors, edited by C Beckers et al., 1989, pp 35-44.

Cai W, Feng D, Fulton R, Siu WC. Generalized linear least squares algorithms for modeling glucose metabolism in the human brain with corrections for vascular effects. Comput Methods Programs Biomed. 2002; 68(1): 1-14. doi: 10.1016/s0169-2607(01)00160-2.

Evans AC. A double integral form of the three-compartmental, four-rate-constant model for faster generation of parametric maps. J Cereb Blood Flow Metab. 1987; 7:S453.

Feng D, Ho D. Parametric imaging algorithms for multicompartmental models - dynamic studies with positron emission tomography. In: Quantification of Brain Function: Tracer Kinetics and Image Analysis in Brain PET, Uemura K, Lassen NA, Jones T, Kanno I (eds.), Elsevier, 1993, pp. 127–137. ISBN: 044489859X.

Feng D, Ho D, Lau KK, Siu WC. GLLS for optimally sampled continuous dynamic system modeling: theory and algorithm. Comput Methods Programs Biomed. 1999; 59(1): 31-43. doi: 10.1016/s0169-2607(98)00099-6.

Gjedde A, Wong DF. Modeling neuroreceptor binding of radioligands in vivo. In: Quantitative imaging, neuroreceptors, neurotransmitters, and enzymes. (Ed. JJ Frost and HN Wagner Jr). Raven Press, New York, 1990, 51-79. ISBN: 978-0881676112.

Gjedde A. Modeling the dopamine system in vivo. In: In vivo imaging of neurotransmitter functions in brain, heart, and tumors. (Ed. Kuhl DE). ACNP Publication 91-2, USA, 1991, 157-179.

Gjedde A. Modelling metabolite and tracer kinetics. In: Molecular Nuclear Medicine: The Challenge of Genomics and Proteomics to Clinical Practice (Eds. Feinendegen LE, Shreeve WW, Eckelman WC, Bahk Y-W, Wagner HN Jr), Springer-Verlag, 2003, pp 121–169. doi: 10.1007/978-3-642-55539-8_7.

Ichise M, Toyama H, Innis RB, Carson RE. Strategies to improve neuroreceptor parameter estimation by linear regression analysis. J Cereb Blood Flow Metab. 2002; 22(10): 1271-1281. doi: 10.1097/01.WCB.0000038000.34930.4E.

Kim SJ, Lee JS, Kim YK, Frost J, Wand G, McCaul ME, Lee DS. Multiple linear analysis methods for the quantification of irreversibly binding radiotracers. J Cereb Blood Flow Metab. 2008; 28(12): 1965-1977. doi: 10.1038/jcbfm.2008.84.

Sederholm K: Using NNLS in multilinear PET problems. TPC modelling reports, 2003. tpcmod0020.pdf.

Thie JA, Smith GT, Hubner KF. Linear least squares compartmental-model-independent parameter identification in PET. IEEE Trans Med Imaging 1997; 16(1): 11-16. doi: 10.1109/42.552051.

Zhou Y, Endres CJ, Brasić JR, Huang SC, Wong DF. Linear regression with spatial constraint to generate parametric images of ligand-receptor dynamic PET studies with a simplified reference tissue model. Neuroimage 2003; 18(4): 975-989. doi: 10.1016/s1053-8119(03)00017-x.

Tags: Compartmental model, NNLS

Updated at: 2023-07-01
Created at: 2023-06-07
Written by: Vesa Oikonen