Analytic solutions of ODEs

The compartmental models for PET data (assuming linear system) can be described in terms of a set of ordinary differential equations (ODEs). The analytic solutions of ODEs can be derived using Laplace transformation, which utilize the operation of convolution, ⊗. Convolution integral of the input function with a response function h(t) gives the concentration in the tissue compartment

$C_{T}(t) = C_{\text{input}}(t) \otimes h(t) = \int_{0}^{t} C_{\text{input}}(\tau) h(t-\tau) d\tau$

The response function is zero when t≤0, and sum of exponentials when t>0; the number of exponentials is the same as the number of tissue compartments.

As an example, for one-tissue compartmental model (1TCM) the solution is:

$C_1(t) = K_1 C_{0}(t) \otimes e^{-k_{2}t}$

The measured PET data is discrete, with relatively long sampling intervals. To achieve unbiased results, convolution operation in practise (linear discrete time convolution method) requires short sample time intervals (Δt). The time-activity curves (discretely sampled data), and the response function, need to be interpolated to even sample times and sampling durations. The value of each sample in calculation of the convolution integral should represent the mean value during the (short) Δt: in practise this means that PET data should be interpolated to the midpoint of the even intervals, and if the response function is an integrable function, then its definite integral during Δt, divided by Δt, should be used. Resulting convolution integral on the other hand should be scaled by Δt; thus Δt is cancelled out. For example, in the case of the 1TCM, the definite integral of the response function between times t₁ and t₂ is:

$\int_{t_1}^{t_2} e^{-k_{2}t} dt = \frac{1}{k_2} \times \left( e^{-k_2 t_1} - e^{-k_2 t_2} \right)$

When convolution integral is calculated at N time points of equal length Δt, that is, at time points ½Δt, Δt+½Δt, 2*Δt+½Δt, ..., (N-1)*Δt+½Δt, the response function values (requiring no further scaling) are calculated as

$\int_{T - \frac{\Delta t}{2}}^{T + \frac{\Delta t}{2}} e^{-k_{2}t} dt = \ \frac{1}{k_2} \times \left( e^{k_2 \Delta t} - 1 \right) e^{-k_2 (T + \frac{\Delta t}{2})}$

Convolution of input TAC with response function consisting of up to three exponential can be performed using convexpf, or specifically applying one-tissue compartment model using conv1tcm.

Alternative solutions that avoid the computationally slow convolution operation include the direct numerical approach (Graham, 1985) and 2nd order Adams-Moulton method based approach (Kuwabara et al., 1993). As an example, the solution of the one-tissue compartmental model ODE is:

$C_{1}(T) = \frac{K_1 \int_{0}^{T} C_{0}(t)dt \ - k_2 \left[ \int_{0}^{T - \Delta t} C_{1}(t)dt + \ \frac{\Delta t}{2} C_{1}(T - \Delta t) \right] }{1 + \frac{\Delta t}{2} k_2}$

Continuing with the example of the one-tissue compartmental model, utilized in the model for [¹⁵O]H₂ PET studies, we simulate renal cortex tissue curve with realistic input function (arterial blood curve) and parameters K₁=3 mL*(min*mL)^-1 and k₂=3.2 min^-1. Arterial data is assumed to be measured from the dynamic image, and therefore have time frames of duration 4 s in the beginning of the scan, and the last time frame is 60 s. Linear discrete time convolution method requires even sample intervals, and to test how the selected interpolation intervals affect the accuracy, three different interpolation intervals are used, 0.1, 1.0, and 4.0 s. For comparison, tissue data is simulated also using the second-order Adams-Moulton method, without any interpolation. Simulations can be done using programs conv1tcm and sim_3tcm, and batch file is available in GitLab. The simulated TACs are shown in Figure 1a, with only the initial part of the data in Figure 1b.

1TCM simulation using convolutions - initial phase of TAC

Figure 1a and 1b. Renal cortex curve is simulated from arterial blood data using the one-tissue compartmental model applying convolution method and Adams-Moulton method. Convolution method is used with three different sample time intervals, 0.1, 1, and 4 seconds. When observing the whole data range (a), all methods seem to lead to similar result. When looking closer to just the initial phase of the simulated tissue curve (b), it is easier to see that with convolution method with long interpolation intervals leads to biased simulation, and with shorter sampling intervals the simulated tissue curve becomes closer to the curve produced by Adams-Moulton used without any interpolation.

Literature

Graham MM. Model simplification: complexity versus reduction. Circulation 1985; 72(5 Pt 2):IV63-IV68. PMID: 4053329.

Kuwabara H, Cumming P, Reith J, Léger G, Diksic M, Evans AC, Gjedde A. Human striatal L-DOPA decarboxylase activity estimated in vivo using 6-[¹⁸F]fluoro-DOPA and positron emission tomography: error analysis and application to normal subjects. J Cereb Blood Flow Metab. 1993; 13: 43-56. doi: 10.1038/jcbfm.1993.7.

Thompson WJ. Computing for Scientists and Engineers - A Workbook of Analysis, Numerics, and Applications. Wiley, 1992. ISBN 0-471-54718-2.

Tags: Modeling, Compartmental model, Simulation, ODE, Convolution

Updated at: 2023-06-14
Created at: 2018-03-10
Written by: Vesa Oikonen

Analytic solutions of ODEs

See also:

Literature