Exponential functions as PET model input

Exponential function has a relationship with compartmental models. In pharmacokinetic one-compartment model drug clearance can be estimated by fitting a single exponential function to the decreasing curve of the drug plasma concentration. Sum of exponential functions can be used to describe the kinetics of models with more compartments.

$\begin{equation} f(t) = \sum_{i=1}^{N} A_i e^{-\lambda_i t} \end{equation}$

In pharmacokinetic analysis the appropriate N can be selected using Akaike Information Criterium (Yamaoka et al., 1978).

Sums of exponential functions (Figure 1) can also be used to fit the decreasing tail of bolus plasma or blood curves collected during a PET study for estimating the plasma pharmacokinetics of the radiotracer.

Examples of exponential functions — **Figure 1.** On the left, examples of three exponential functions, with parameters *A₁=5* and *λ₁=0.2* (red), *A₂=0.5* and *λ₂=0.01* (blue), and *A₃=0.2* and *λ₃=0.0001* (green). On the right, the sum of these three exponential functions.
In this example all time constants (λs, eigenvalues) are negative, and therefore decreasing with time. At time zero *f(0)=A_i*. If *λ=0*, then function gets value *A_i* at all times, otherwise function values approach zero. The third exponential has very low λ, and therefore the sum function stays at ∼*A₃=0.2* during this time range.

Example of sum of exponential functions — **Figure 1.** On the left, examples of three exponential functions, with parameters *A₁=5* and *λ₁=0.2* (red), *A₂=0.5* and *λ₂=0.01* (blue), and *A₃=0.2* and *λ₃=0.0001* (green). On the right, the sum of these three exponential functions.
In this example all time constants (λs, eigenvalues) are negative, and therefore decreasing with time. At time zero *f(0)=A_i*. If *λ=0*, then function gets value *A_i* at all times, otherwise function values approach zero. The third exponential has very low λ, and therefore the sum function stays at ∼*A₃=0.2* during this time range.

Monoexponential extrapolation of radiopharmaceutical clearance curves is a common procedure (Lassen & Sejrsen, 1971). For extrapolating the input curve, extrapol can be used. However, better approach usually is to write a script which first fits 2 or 3 exponential function with fit_exp to the input data starting at some time after the peak is passed, or automatically determined number of exponentials after peak with fit_dexp; then computes the fitted data curve at desired time points using fit2dat, and finally combines the original data and fitted tail using taccat.

Arterial curves after bolus injection can be very concave, and decent fitting requires N≥3. Instead of that, residence times can be approximated by gamma distribution, which can be interpreted in terms of recirculatory pharmacokinetic model (Weiss, 1983), or with a form of Weibull distribution (Piotrovskii, 1987a, 1987b), using either of the functions

$\begin{equation} f(t) = A \: t^{-\beta} \: e^{-\lambda t} \end{equation}$ $\begin{equation} f(t) = A \: \beta \: \lambda \: t^{-(\beta - 1)} \: e^{-\lambda t^{\beta}} \end{equation}$

instead of one of the exponential terms. In these functions 0<β<1. By setting β to zero or one, respectively, the functions become normal exponentials.

Actual arterial PET input curves, following intravenous bolus administration of radiotracer, typically start with a period of zero concentration because of the delay that it takes for the bolus to reach the sampling site. After that follows a steep rise, then the concentration curve passes through a maximum, and then undergo a slower decay. While the descending part of the bolus curve can be fitted using the sum of decaying exponentials, the ascending phase (before peak time, T_peak) can be fitted using a line (Parsey et al., 2005; Vriens et al., 2009):

$\begin{equation} f(t) = \left\{ \begin{array}{l l} 0 & \quad \text{if } t < -\frac{b}{a}\\ a \times t + b & \quad \text{if } -\frac{b}{a} \leq t < T_{peak}\\ \sum_{i=1}^{3} {A_i} e^{-\lambda_i (t-T_{peak})} & \quad \text{if } t \geq T_{peak} \end{array} \right. \end{equation}$

, or, assuming that ascending phase starts at time 0 (Su et al., 1994):

$\begin{equation} f(t) = \left\{ \begin{array}{l l} (\sum_{i=1}^{N} {A_i}) \frac{t}{T_{peak}} & \quad \small{\text{if } t \leq T_{peak} }\\ \sum_{i=1}^{N} {A_i} e^{-\lambda_i (t-T_{peak})} & \quad \small{\text{if } t \geq T_{peak} } \end{array} \right. \end{equation}$

Ascending phase has also been fitted using a 2nd or 3rd order polynomial function (Eberl et al., 1997; Bentourkia et al., 1999), which however can lead to negative function values. Non-negativity can be assured by fitting the ascending phase using instead function (Wong et al., 2006):

$f(t) = a \: (1 - \frac{1}{1 + b \: t^3})$

Based on a compartmentalized model of radiopharmaceutical behaviour in the circulatory system, Feng et al. (1993a and 1993b) have proposed a formulation of four exponential functions for fitting PET radiotracer PTACs which include both ascending and descending phase in the measured data, and the initial zero phase before the radiopharmaceutical appearance time, T_ap:

$\begin{equation} f(t) = \left\{ \begin{array}{l l} 0 & \quad \small{\text{if } t \leq T_{ap} } \\ A_1 e^{-\lambda_1 (t-T_{ap})} + A_2 e^{-\lambda_2 (t-T_{ap})} + A_3 e^{-\lambda_3 (t-T_{ap})} & \quad \small{\text{if } t > T_{ap}} \\ - (A_1 + A_2 + A_3) e^{-\lambda_4 (t-T_{ap})} \end{array} \right. \end{equation}$

Eigenvalues (λ₁-λ₄) are ≥0. Examples of curves based on that function are shown in Figure 2. Functions derivative is

$\begin{equation} f^{\prime}(t) = \left\{ \begin{array}{l l} 0 & \quad \small{\text{if } t \leq T_{ap} } \\ - A_1 \lambda_1 e^{-\lambda_1 (t-T_{ap})} - A_2 \lambda_2 e^{-\lambda_2 (t-T_{ap})} - A_3 \lambda_3 e^{-\lambda_3 (t-T_{ap})} & \quad \small{\text{if } t > T_{ap} } \\ + (A_1 + A_2 + A_3) \lambda_4 e^{-\lambda_4 (t-T_{ap})} \end{array} \right. \end{equation}$

To ensure that the function does not give negative (non-physiological) values at the start, that is, f'(T_ap)>0, this condition has to be met:

$(A_1 + A_2 + A_3) \: \lambda_4 - A_1 \lambda_1 - A_2 \lambda_2 - A_3 \lambda_3 \gt 0$

Examples of sum of exponential functions for bolus studies — **Figure 2.** Examples of curves that can be produced with the four-exponential function.
The parameters of the curves are *A₁=0.1*, *λ₁=0.0001*, *A₂=3*, *λ₂=0.01*, *A₃=1*, *λ₃=0.1*, and *λ₄=0.5* (red); *A₁=1*, *λ₁=0.001*, *A₂=2*, *λ₂=0.2*, *A₃=5*, *λ₃=0.02*, and *λ₄=0.1* (blue); *A₁=4*, *λ₁=0.01*, *A₂=4*, *λ₂=0.1*, *A₃=2*, *λ₃=0.001*, and *λ₄=0.05* (green). For all curves *T_ap=50*.

Two of the eigenvalues can be paired, leading to this function.

Fitting a multiexponential function is an ill-posed problem and not that easy to solve as it might first appear (Istratov & Vyvenko, 1999). For fitting the descending exponential part, initial estimates have traditionally been estimated using curve stripping techniques (Dunne, 1986); Initial parameters for unconstrained multiexponential functions can be calculated using linear methods as suggested by (Jean Jacquelin); this method can be tested with program llsqe3. spectral analysis may be a more robust technique for estimating the parameters of the exponential functions (Bentourkia et al., 1999).

Especially in DCE-MRI studies, but also in PET field, "biexponential" may refer to the sum of two exponentials (McGrath et al., 2009), a function presented above with N=2, or, "biexponential AIF model", which is actually a compartmental model used to calculate realistic BTAC from a rectangular (boxcar) function that represents the infusion of the contrast agent or radiotracer.

Infusion

To simulate the radiotracer infusion, of duration T_in, the boxcar function can be convolved with the sum of exponentials function (Eq 1). TAC can be convolved with sum of exponentials using convexpf. Convolution operation can here be replaced by subtracting integrals of the exponentials. The integral of decaying exponential function (λ_i > 0) is

$\begin{equation} g(t) = \int_0^T{f(u)du} = \sum_{i=1}^{N} \frac{A_i}{\lambda_i} (1 - e^{-\lambda_i T}) \end{equation}$

and the subtraction is calculated as

$\begin{equation} f(t) = g(t) - g(t - T_{in}) \end{equation}$

The function can be calculated in three parts:

$\begin{equation} f(t) = \left\{ \begin{array}{l l} 0 & \quad \small{\text{if } t \leq T_{ap} } \\ \sum_{i=1}^{N} \frac{A_i}{\lambda_i} (1 - e^{-\lambda_i (t-T_{ap})}) & \quad \small{\text{if } T_{ap} < t < (T_{ap} + T_{in}) } \\ \sum_{i=1}^{N} \frac{A_i}{\lambda_i} (e^{-\lambda_i (t-T_{ap}-T_{in})} - e^{-\lambda_i (t-T_{ap})}) & \quad \small{\text{if } t \geq (T_{ap} + T_{in}) } \end{array} \right. \end{equation}$

In pharmacokinetics and simulations the plasma concentration during and after intravenous infusion have been described using this (or with notation differences) function (Wagner, 1976; Madsen et al., 1993; Oakes et al., 1999). This function, optionally with a few extensions, can be fitted with program fit_sinf. T_ap is the appearance time of radioactivity in the blood or plasma, and T_in is the infusion time. N should usually be set to 2 or 3, but in some situations higher number of parameters is required.

The integrals of the same exponential functions that were shown in Figure 1 are shown in Figure 3, together with the infusion model calculated with T_ap=0 and T_in=100. As can be seen from the plots, the function models well the initial delay and dispersion that happens when the bolus travels from IV injection site through the heart cavities and lungs to the arterial system, and the sum function can basically handle the recirculation and redistribution. However, the measured arterial curves are further dispersed in the circulation, and unless curve is dispersion corrected, it cannot be well fitted with this function. If dispersion time constant (τ) is known, then the effect of dispersion can be added to the function in fit_sinf; this helps to get markedly better fits, and allows noise-free dispersion correction.

Examples of integrals of exponential functions — **Figure 3.** On the left, the integrals of three exponential functions (red, blue, and green), and their sum (purple). Exponential functions are the same as in Figure 1. On the right, the plasma curves simulating radiotracer infusion of length *T_in=100*, otherwise again calculated using the same exponential functions.

Example of infusion curves from sum of exponential functions — **Figure 3.** On the left, the integrals of three exponential functions (red, blue, and green), and their sum (purple). Exponential functions are the same as in Figure 1. On the right, the plasma curves simulating radiotracer infusion of length *T_in=100*, otherwise again calculated using the same exponential functions.

Two exponentials are sufficient to fit adequately the blood TAC from RV cavity in myocardial [¹⁵O]H₂O studies. Program fit_wrlv can be used to fit the two-exponential function to RV and LV curves simultaneously, fitting also the delay time and dispersion constant between the RV and LV blood TACs.

Recirculation

The integral of exponential function has been often used for describing the recirculation phase of the plasma TACs in pharmacokinetics and also in PET. The plasma curve is fitted to a function that is the sum of exponential function and its integral

$\begin{equation} f(t) = \sum_{i=1}^{N} \left( A_i e^{-\lambda_i t} + c \times \frac{A_i}{\lambda_i} (1 - e^{-\lambda_i t}) \right) \end{equation}$

, where c is a scaling factor. The exponential functions that were shown in Figure 1 and their scaled integrals (Figure 3, left side) are summed and shown in Figure 4. Similar and more versatile curve shapes can be achieved without the integral function if more exponentials are added, but that increases also the number of function parameters which complicates the fitting. The same approach can be applied to any function.

Examples of exponential functions plus their scaled integrals — **Figure 4.** Examples of curves that can be produced by summing the exponential functions shown in Figure 1 and their scaled integral functions shown in the left panel of Figure 3. The scaling factor *c=0.001*, except for the third function (green) it is set to zero to prevent getting a constantly increasing curve. The purple curve represents the sum of the three functions.

References

Feng D, Huang S-C, Wang X. Models for computer simulation studies of input functions for tracer kinetic modeling with positron emission tomography. Int J Biomed Comput. 1993; 32: 95-110. doi: 10.1016/0020-7101(93)90049-C.

Feng D, Wang Z. A three-stage parameter estimation algorithm for tracer concentration kinetic modelling with positron emission tomography. Proceedings, 1991 American Control Conference, vol 2 (1991): 1404-1405. doi: 10.23919/ACC.1991.4791609.

Lüdemann L, Sreenivasa G, Michel R, Rosner C, Plotkin M, Felix R, Wust P, Amthauer H. Corrections of arterial input function for dynamic H₂¹⁵O PET to assess perfusion of pelvic tumours: arterial blood sampling versus image extraction. Phys Med Biol. 2006; 51: 2883-2900. doi: 10.1088/0031-9155/51/11/014.

Wagner JG. Linear pharmacokinetic equations allowing direct calculation of many needed pharmacokinetic parameters from the coefficients and exponents of polyexponential equations which have been fitted to the data. J Pharmacokin Biopharm. 1976; 4(5): 443-467. doi: 10.1007/BF01062831.

Appendix

Definite integrals of exponential functions:

$\begin{equation} \int_{0}^{T}e^{-k t} dt = \frac{1 - e^{-k T}}{k} \nonumber \end{equation}$

$\begin{equation} \int_{0}^{T}(1 - e^{-k t}) dt = T - \frac{1 - e^{-k T}}{k} \nonumber \end{equation}$

Tags: Input function, Fitting, Extrapolation, Exponential function, Biexponential

Updated at: 2021-12-10
Created at: 2016-08-08
Written by: Vesa Oikonen

Exponential functions as PET model input

Infusion

Recirculation

See also:

References

Appendix