Fourier domain closed-form formulas for estimation of kinetic parameters in reversible multi-compartment models.

BACKGROUND: Compared with static imaging, dynamic emission computed tomographic imaging with compartment modeling can quantify in vivo physiologic processes, providing useful information about molecular disease processes. Dynamic imaging involves estimation of kinetic rate parameters. For multi-compartment models, kinetic parameter estimation can be computationally demanding and problematic with local minima.
METHODS: This paper offers a new perspective to the compartment model fitting problem where Fourier linear system theory is applied to derive closed-form formulas for estimating kinetic parameters for the two-compartment model. The proposed Fourier domain estimation method provides a unique solution, and offers very different noise response as compared to traditional non-linear chi-squared minimization techniques.
RESULTS: The unique feature of the proposed Fourier domain method is that only low frequency components are used for kinetic parameter estimation, where the DC (i.e., the zero frequency) component in the data is treated as the most important information, and high frequency components that tend to be corrupted by statistical noise are discarded. Computer simulations show that the proposed method is robust without having to specify the initial condition. The resultant solution can be fine tuned using the traditional iterative method.
CONCLUSIONS: The proposed Fourier-domain estimation method has closed-form formulas. The proposed Fourier-domain curve-fitting method does not require an initial condition, it minimizes a quadratic objective function and a closed-form solution can be obtained. The noise is easier to control, simply by discarding the high frequency components, and emphasizing the DC component.

Introduction

Dynamic emission computed tomographic imaging can measure the kinetics of the tracer’s distribution and exchange with body tissues. Using quantitative analysis techniques such as compartment modeling [1-4], dynamic imaging can quantify in vivo physiologic and metabolic processes, providing more information regarding underlying molecular disease processes than can be obtained from static imaging. However, the fitting of compartment models to dynamic imaging data can be computationally demanding and can have problems with local minima.

A number of techniques for kinetic parameter estimation have been studied and are in use today, generally offering a tradeoff between computation time, robustness of fit, and flexibility with differing sets of assumptions. Perhaps the most robust—but also most computationally demanding—approach for estimating individual rate parameters for multi-compartment models is classic nonlinear least-squares estimation [5]. Here, a least-squares or weighted least-squares objective function is iteratively minimized to obtain best-fit kinetic parameter estimates. Numerous curve-fitting algorithms have been investigated for this application, including derivative-based or downhill simplex methods [5], ridge-regression [6-11], simulated annealing [12], and even discretized exhaustive search paradigms. The best results are obtained when accurate weighting is applied, though determination of the best weights is itself a challenging problem that depends on many variables including the reconstruction algorithm used [13]. Unfortunately, the least-squares objective function for multi-compartment models with noisy data is ill-formed, containing broad shallow valleys and (potentially) local minima. As such, careful implementation of the curve-fitting algorithm with extensive iteration and handling of local minima is required to confidently find the global minimum. This results in computational-demands that, while reasonable for fitting individual time-activity curves, may become impractical for voxelwise parametric imaging where millions of fits need to be performed.

The conventional least-squares curve-fitting method is to match the measured time-activity-curve with the calculated solution of the ordinary differential equations. The solution is in the form of a non-linear function of time and the kinetic parameters. It is neither in the ODE form or the integral form. Much of the difficulty in the curve-fitting approaches is due to the presence of nonlinear terms in the solution to the compartment modeling equations. Significant efforts have been made to “linearize” the problem, with milestones including the weighted integration method [14-16], linear least-squares (LLS) and generalized linear least-squares (GLLS) methods [17-20]. Such methods are based on integrating the compartment modeling equations to obtain linear systems of equations relating the rate parameters (or combinations thereof) to integrals of the time-activity curves and blood input functions. Though the image noise is not correlated between timeframes, integrating the time-activity curves in time gives equation errors that are not statistically independent, and that can bias the results (e.g. LLS); this correlation-induced bias can be reduced or eliminated by the iterative auto-regressive filtering technique of GLLS. More recently, approaches employing temporal basis functions to reduce noise and improve parameter estimation have also been studied [21-26]. Performance comparisons of these methods have been performed by Feng et al.[27] and more recently by Dai et al.[28]. Overall, nonlinear least-squares was found to provide the most robust parameter estimates, though this approach is also the most computationally intensive and initial condition dependent. Of the fast approaches, GLLS performed well for lower noise data, but may exhibit large bias and poor precisions when the noise level is high. Certain basis function approaches appear promising, though they are currently slower than GLLS and less thoroughly investigated.

This work offers a new perspective to the compartment model fitting problem where Fourier linear system theory is applied to derive closed-form formulas for estimating kinetic parameters for the two-tissue compartment model. This paper is an extension of our conference presentation at 2011 IEEE Medical Imaging Conference [29]. It builds on our related work developing closed-form solutions in the time domain [30], where the continuous-time differential equation is transformed into a discrete-time difference equation by using the exact continuous-time expression of the compartment activity. This time-domain approach is fundamentally different from time-domain curve-fitting approaches, and thus will likely have very different characteristics in terms of accuracy, speed, noise propagation and precision.

This paper introduces Fourier-domain closed-form solutions for the estimation of kinetic parameters for the two-tissue compartment model with 4 rate parameters. The similarity of the time-domain method and this Fourier-domain method are in the closed-form strategy to find the optimal solution of an objective function. Implementation procedures of the approach are discussed, and initial computer simulations are performed demonstrating application of the technique. The unique feature of the proposed method is that the formulation is in the Fourier domain, where the high-frequency noisy components can be readily discarded, and exponential functions never appear.

Methods and results

This section presents the standard method of non-linear fitting to estimate the kinetic parameters, points out its potential deficiencies of the dependency on the initial conditions, and develops a closed-form Fourier domain method that is independent of the initial conditions. The result of the proposed Fourier domain estimation method is close to the true solution and can be used as the initial condition for the standard iterative method for refining.

The standard iterative curve-fitting method

Figure 1 shows a generic two-tissue compartment model with four rate constants: K1, k2, k3 and k4. The model can be described by two first-order differential equations:

Figure 1

A general two-compartment-model.

The above two equations are referred to as the state equations in system theory. The time-activity curve for the blood input is B(t). There are two tissue compartments C1(t) and C2(t); these two compartments are not individually accessible and their time-activity curves cannot be individually measured; however, the sum of them, C(t) defined as

can be measured. Equation (15) is called the output equation in system theory. Combining (1), (2) and (3) yields a second order differential equation:

with

The solution of (4) can be expressed as two convolution terms:

where

The standard iterative parameter estimation method is to minimize the objective function:

which measures the discrepancy between the measured tissue time-activity-curve C(t) and the estimated time-activity-curve C(t) evaluated by (6).

More precisely, , where C(t) is the measured activity count at time t and w is the associated reciprocal variance. Also note C(t) need not be in the form of a convolution integral. We can alternately use an ODE-solver in conjunction with curve-fitting to compute C(t) directly from the (1), (2), and (3) as needed in each curve-fitting iteration.

Since the gradient of the objective function F in (11) is non-linear with respect to the unknown parameters: K1, k2, k3 and k4, the result of the iterative algorithm will in general depend on the initial conditions. This phenomenon will be illustrated by two examples in the next section. After the examples, we will present a Fourier domain estimation method to solve the problem of dependency on the initial condition. The Fourier domain estimation method is the main contribution of this paper.

Numerical examples of the standard iterative curve-fitting method

In the following computer simulations, the blood input functions were in the form of [31]:

where A1 = 22.24, A2 = 8.36, A3 = 0.10, λ1 = 21.28 min-1, λ2 = 7.71 min-1, and λ3 = 0.37 min-1. The blood and tissue time activity curves were non-uniformly sampled as 30×5.sec., 20×10 sec., 10×30 sec., 10×60 sec., 10×150 sec., and 3×300 sec. At each sampling time, the activity was integrated over a 60-second time window. The total scanning time was about one hour (i.e., 3650 seconds).

Measurements as “30×5 sec” means that 30 samples are taken and the samples are 5 sec apart between adjacent samples. Each sample is a continuous time integral of the time-activity-curve, and the integration time interval is 60 sec. The upper limit of the time interval is the sampling time instant. Therefore, the sampled signal is a collection of overlapped time integrals of the actual time-activity-curve. Even when the gap between adjacent samples is 300 sec, the integration time interval is still 60 sec, and in this case there is no overlap for time integral intervals for the samples.

The purpose of time integration is to reduce noise. Modern nuclear medicine scanners can acquire data in list-mode, in which data are acquired continuously with a non-discretized time stamp for each event. After the list-mode data are acquired, the discrete samples are formed for image reconstruction. When the discrete samples are formed, the time integral (i.e., count summation) is performed.

In the numerical examples, we chose K1 = 0.4 min-1, k2 = 0.3 min-1, k3 = 0.2 min-1 and k4 = 0.1 min-1. The k values used in the computer simulations were tailored from [27] and [32], where a wide range of the k values are reported for different patient studies. For example, K1 can be anywhere from 0.1 to1.2, and k2 can be from 0.1 to 0.8, and so on. The values of k3 and k4 are smaller. Scaled Gaussian noise, N(0,1) (mean = 0, standard deviation = 1), was added to the noiseless data C(t). The noise scaling factor was

where T1/2 (=110 min) was the half-life of the isotope, and the proportional constant α = 0 corresponds to the noiseless case. This noise model is suggested in [31]. Typical noisy time activity curves C(t) are shown in Figure 2 for these three α values. In reality, the images are reconstructed at the pre-selected non-uniformly distributed time instances, and the non-uniformly distributed blood and tissue time-activity-curves are obtained.

Figure 2

Time-integrated time-activity curves () with 3 different noise levels. The total integration interval is one minute.

The solution of (1)-(3) corresponding to is denoted as Ĉ(t). Let B(t) and C(t) denote the corresponding curves after time integration over the 60 sec interval. {, Ĉ(t)} and {B(t), C(t)} satisfy the second order differential equation (4) based on the theory to be discussed in Data integration effects section.

The standard least-squares fitting method was applied to the model (6) using MATLAB® and used the built-in function “nlinfit.” The MATLAB’s built-in function nlinfit allows negative outcomes. MATLAB also has another built-in function lsqnonlin that enforces the non-negativity constraint on the outcomes and is widely used for non-negative parameter fits.

Some numerical examples are shown in Tables 1, 2, and indicate the importance and sensitivity of the iterative parameter estimation algorithm in selection of a proper initial condition. In Table 1, only one noise realization was used for each case, while in Table 2, 250 noise realizations were used for each case. If a bad initial condition is chosen, the results could be totally wrong even under the ideal noiseless situation. In order to overcome this problem, a closed-form Fourier domain estimation is developed in the next section. In Table 1, the true k parameters are not the same as the k parameters estimated with noiseless (α = 0) data. In the noiseless (α = 0) case, both B(t) and C(t) are error free. If the initial conditions are not properly provided, the estimated k parameters from the noiseless data are still far away from the true parameters.

Table 1

The results of a standard iterative curve-fitting method with a good and a bad initial condition (using 1 noise realization for each case)

Case	Bad initial condition				Good initial condition
	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)
Initial	0.70	0.70	0.70	0.70	0.45	0.35	0.25	0.15
α = 0	0.14	−1.62	−5.10	−2.10	0.40	0.30	0.20	0.10
α = 0.1	0.13	−1.65	−5.14	−2.090	0.4	0.28	0.18	0.10
α = 0.4	0.12	−1.64	−5.03	−2.038	0.391	0.26	0.15	0.09
True	0.40	0.30	0.20	0.10	0.40	0.30	0.20	0.10

Table 2

The results of a standard iterative curve-fitting method with a good and a bad initial condition (using 250 noise realizations for each case)

Case	Bad initial condition				Good initial condition
	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)
Initial	0.70	0.70	0.70	0.70	0.45	0.35	0.25	0.15
α = 0	0.14	−1.62	−5.10	−2.10	0.40	0.30	0.20	0.10
α = 0.1	0.15 ± 0.01	−1.61 ± 0.07	−5.08 ± 0.16	−2.10 ± 0.05	040 ± 0.00	0.30 ± 0.01	0.20 ± 0.02	0.10 ± 0.00
α = 0.4	0.14 ± 0.05	−1.62 ± 0.34	−5.11 ± 0.78	−2.10 ± 0.27	040 ± 0.02	0.31 ± 0.07	0.21 ± 0.06	0.10 ± 0.01
True	0.40	0.30	0.20	0.10	0.40	0.30	0.20	0.10

The Fourier domain estimation method

By taking the Fourier transform of (4), we have

where and are the Fourier transform of C(t) and B(t), respectively, and ω is the frequency variable. Let ω = 0, (14) immediately gives a DC (direct current) gain relationship:

In the time domain, (15) is expressed as

If the noise in B(t) and C(t) has a zero mean value over time, (16) is not affected by noise.

We now introduce some short-hand notations:

Using these short-hand notations, Eq. (14) becomes

Substituting (15), namely,, into (18) eliminates x4, and an objective function H can be formulated as:

where

In order to find the minimum of the objective function H, we set the partial derivatives of H to 0 and obtain a set of linear equations as follows.

where

In (22) and (24), “Re” denotes the operation of “taking the real part,” and the inner-product is defined by

A closed-form solution can readily be obtained as

Finally, the two-compartment model kinetic parameters are estimated as

TAC extrapolation

If a prolonged data acquisition is not practical, the un-measured “tail” activity curve must be estimated by using extrapolation methods. Since the tail decays monotonically and does not contain any peaks, the tail estimation error is well controlled.

The proposed method depends heavily on the DC gain (or the VD value). If the measurement of the TAC C(t) is terminated before it reaches zero, the VD value still can be estimated [33]. An alternative approach is to estimate the unmeasured C(t) decay trend as follows.

Here, we suggest one data extrapolation method that uses a summation of a family of exponential functions with different decay constants to approximate the truncated “tail,” under the constraint that the activity curves tend to zero when time approaches infinity. The truncated TAC C(t) can be approximated by

where t is the time when data acquisition stops and λ (n = 1, 2, …, 10) are a family of decay constants and are 0.01, 0.02, …, 0.10 in our illustration examples. The motivation to use a family of exponential decays to approximate the unmeasured data is to avoid any bias towards any particular decay rate which is unknown. In computer simulations, C(t) was the average value of the last 3 measurements.

Implementation

The computer procedure is discussed in this section. For the two-compartment model, the steps are:

Step 1: Acquire the blood input function B(t) and the compartment time-activity curve C(t), for n = 0, 1, 2, …, N-1. Here the sampling interval can be non-uniform and N is the total number of samples.

Step 2: Extrapolate the time-activity-curve using the method described in TAC extrapolation section. Linear interpolate the data so that the resultant data sets have a constant sampling interval T. The total number of the new data points becomes a much larger number N.

Step 3: Evaluate and .

Step 4: Use a DFT (Discrete Fourier Transform) or FFT (Fast Fourier Transform) computer routine to calculate the DFT functions and .

Step 5: Calculate , , , and .

Step 6: Calculate the inner products that appear in (22) and (24) using the inner product definition (25) and a user selected cutoff index k , corresponding to the continuous cutoff frequency ω.

Step 7: Form the matrices and as shown in (22) and (24).

Step 8: Solve for using (26).

Step 9: Finally, use (27) to obtain estimated kinetic parameters K1 ~ k4.

The proposed Fourier domain estimation method is in closed-form and does not need an initial condition.

We now present the results of our algorithm applied to the same computer generated data as in the standard iterative curve-fitting method, the outputs of the proposed Fourier domain are shown in Table 3 as the “Stage 1”. In Table 3, only one noise realization is used for each case, while in Table 4, 250 noise realizations are used for each case. Ten low frequency components were used in the Fourier domain calculation, that is, the cutoff index k was chosen to be 10. The non-uniformly sampled, one-hour time-activity-curves were extrapolated into two-hour curves, and then linearly interpolated with a constant re-sampling interval of 5 seconds.

Table 3

The results of the combined method

Case	Stage 1: Fourier method output				Stage 2: Final iterative method output
	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)
α = 0	0.39	0.27	0.17	0.09	0.40	0.30	0.20	0.10
α = 0.1	0.39	0.27	0.11	0.09	0.40	0.28	0.18	0.10
α = 0.4	0.38	0.23	0.14	0.09	0.42	0.35	0.28	0.12
True	0.40	0.30	0.20	0.10	0.40	0.30	0.20	0.10

The output of the proposed closed-form method is used as the initial condition of the iterative method. One noise realization is used.

Table 4

The results of the combined method

Case	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4 (min-1)	K1(k3 + k4) /(k2k4)	K1k3/(k2k4)	k3/k4
α = 0	0.40	0.30	0.20	0.10	4.00	2.67	2.00
α = 0.1	0.40 ± 0.00	0.30 ± 0.01	0.20 ± 0.01	0.10 ± 0.00	4.00 ± 0.02	2.67 ± 0.05	2.01 ± 0.11
α = 0.4	0.40 ± 0.02	0.31 ± 0.07	0.21 ± 0.06	0.10 ± 0.01	4.01 ± 0.08	2.68 ± 0.19	2.07 ± 0.46
True	0.40	0.30	0.20	0.10	4.00	2.67	2.00

The output of the proposed closed-form method is used as the initial condition of the iterative method. The number of noise realization is 250.

Some macro-parameters have special clinical significance. For example, VD indicates volume-of-distribution, (K1k3)/(k2 + k3) reflects the net uptake, and (K1k3)/( k2k4) or k3/ k4 are related to the binding potential. Therefore, the estimation of some macro-parameters is also reported in the results in Table 4.

It is noticed from Stage 1 of Table 3 that the results are slightly biased even in the noiseless situation. This is because the data are discretely sampled and the continuous Fourier integration must be approximated by a discrete Fourier transform. Another source of error comes from data extrapolation which is required in this Fourier-domain closed-form method. The requirement is that the Fourier integration interval must contain all non-zero activities. Therefore, it is expected that the traditional iterative non-linear fitting method should have less bias, because neither data extrapolation nor interpolation are necessary.

We have mentioned that when a proper initial condition is provided, the standard least-squares fitting usually gives satisfactory results. However, if an inappropriate initial condition is provided, the algorithm may converge to a wrong solution. When the closed-form method was used to obtain a rough estimation and the iterative method was used to fine tune the solution, much better results were obtained as shown in Table 3.

Discussion

Some issues related to the proposed closed-form formula are discussed below.

Data integration effects

Averaging of the time-activity curve (TAC) over a time frame does not affect the accuracy of the equations. Because the TAC C(t) and the input function B(t) are related by a linear differential equation, for example (4), C(t) can be expressed as a convolution with B(t):

where the kernel H(t) is determined by the system of differential equation. Averaging C(t) over a time frame is to convolve C(t) with a box-car window function or unit-step window function W(t):

Combining the above two convolution expressions yields

which implies that if we integrate the TAC C(t) and the input function B(t) over the same time interval, the integrated functions still satisfy the same differential equations.

Selection of the cut-off frequency

The selection of the cut-off s based on trial-and-error, and is dependent on noise level. Ideally speaking, the cut-off frequency should be chosen as high as possible to include as much information. However, the high frequency components contain more noise than the low frequency components. If the cut-off frequency is selected too high, the estimation will be affected by the noise fluctuation, resulting in different solutions with different noise realizations. If the cut-off frequency is selected too low, even though the estimation is not sensitive to noise, it may contain large bias because the system is not adequately represented. The effect of the cut-off frequency selection is illustrated by the examples in Table 5, where a two-compartment model is considered with noise level α = 0.4.

Table 5

Estimation results for the two-compartment model parameters , , , and with noise level α = 0.4, using different cut-off frequencies and 250 noise realizations

Cut-off frequency index kc	K1 (min-1)	k2 (min-1)	k3 (min-1)	k4(min-1)	K1(k3 + k4) /(k2k4)	K1k3/(k2k4)	k3/k4
8	0.36 ± 0.05	0.22 ± 0.10	0.11 ± 0.09	0.07 ± 0.03	4.08 ± 0.07	2.24 ± 0.55	1.43 ± 0.77
10	0.37 ± 0.05	0.24 ± 0.11	0.13 ± 0.08	0.08 ± 0.03	4.07 ± 0.07	2.35 ± 0.49	1.58 ± 0.77
15	0.37 ± 0.04	0.23 ± 0.09	0.11 ± 0.08	0.07 ± 0.03	4.07 ± 0.07	2.27 ± 0.51	1.45 ± 0.71
20	0.37 ± 0.03	0.22 ± 0.08	0.10 ± 0.08	0.06 ± 0.03	4.08 ± 0.07	2.24 ± 0.51	1.39 ± 0.65
True	0.40	0.30	0.20	0.10	4.00	2.67	2.00

Extension to one-compartment or multi-compartment models

In general, if we have N, where N can be 1, compartments for the tissue model, we have a system of N first-order differential equations (which are called state equations) to describe the kinetics [34]. Let C be a vector that contains all N compartments, A be an N×N matrix, D be an N×1 matrix and E be a 1×N matrix. The system of differential equations can be expressed in a matrix form

The measurable activity is described by the output equation:

where the C(t) on the left-hand-side is a scalar. Using the Fourier transform, the system’s transfer function can be obtained as [35]

namely,

Let ω = 0, (35) immediately gives a DC gain relation:

Similar to the case in the Methods and Results section, the coefficients x1 … x are related to the kinetic parameters. The parameter estimation problem can be achieved by minimizing an objective function similar to (19). By using the DC gain expressed in (36), the total number of unknowns is 2 N-1. The unknown coefficients are obtained by solving a least-squares problem to approximate (35). Finally the kinetic parameters are estimated using relations between the coefficients x1 … x and the kinetic parameters. The proposed method becomes less effective as N increases, partially due to the non-linear relationship between the k parameters and x parameters. For a large N, there may not be a closed-form expressions for the k parameter for a set of x’s.

Consideration of the input function contamination effect

In the above discussion, we assume that the quantity C(t) can be measured. In reality, the measured C(t) may be contaminated by the input function B(t). The contamination of the tissue by the blood input function is considered in this section. The contamination of the blood by the tissue is more problematic because it may result in a non-identifiable problem. The system’s transfer function or the impulse response function should be revised to reflect this effect. The revised impulse response function becomes

where δ(t) is the Dirac delta function, and f is the fraction of the contamination. Using the two-compartment model as an example, the revised transfer function is

If we follow the same steps as outlined in Methods and Results section and treat f as another unknown parameter, closed-form formulas for the kinetic parameters as well as f can be obtained.

Relation to other linear methods

Generally speaking, a frequency-domain method is equivalent to a time-domain method. Using a cut-off in the frequency-domain is equivalent to performing smoothing in the time-domain. However, the discretization can introduce some complications. The time domain approximation of the derivative B’(t) by finite difference is not quite equivalent to in the Fourier domain. All linear methods such as LLS and GLLS have their different ways to converting differential equations into algebraic equations. The goal is to remove the (continuous) derivative operator, which cannot be implemented directly with sampled data. In the LLS and GLLS method, the derivative operator is removed by integrating both sides of the differential equation. In the proposed Fourier domain method, the derivative operator is indirectly implemented in the Fourier domain as a multiplication of ikΩ. These two approaches may not be equivalent.

Special treatment of the DC component

We give the DC component special attention because of noise considerations. When the noise is zero-mean, the DC-component is almost un-affected; however, all other frequency components are affected. In a system of equations, we trust the DC relationship the most and enforce it to be a constraint. We discard the noise-heavily-affected components (i.e., the high frequency components). We then weigh the remaining frequency components evenly, as an un-weighted least-square objective function. We control the noise by dividing the frequency components in three categories: Trustworthy (i.e., the DC), somewhat trustworthy (i.e., the low-frequencies), and not-trustworthy (i.e., the high frequencies). In cases where the initial estimate of VD is more problematic than in our current simulations, it is possible to remove the VD constraint and use one more unknown variable in the formulation of the closed-form solution.

We have made another version of the analytical method, in which we do not use the DC gain as a constraint. The DC term is treated almost the same as other low-frequency components but with 1000 times heavier weights. In our simulations, the resultant algorithm was not as stable as the proposed algorithm, in which the DC gain is used as a constraint.

Conclusions

Time-domain curve-fitting is the current state-of-the-art in nuclear medicine kinetic estimation. Due to the non-linear exponential functions, this curve-fitting is sensitive to initial conditions (i.e., the initial solutions) for multi-compartment model parameter estimation problems. The initial condition may make the algorithm converge to a wrong solution. In this paper, a Fourier-domain kinetic estimation method is proposed. The proposed Fourier-domain curve-fitting method does not require an initial condition, it minimizes a quadratic objective function and a closed-form solution can be obtained. The noise is easier to control, simply by discarding the high frequency components, and emphasizing the DC component. The proposed Fourier-domain estimation method has closed-form formulas.

The unique strategy in this paper is to formulate the objective function in the Fourier domain. This strategy has two advantages: The model is linear in terms of transformed variables (18) and it is easier to control noise by discarding high frequency components. The Frequency domain objective function is a quadratic function. To minimize it, the gradients are set to zero. This results in a system of linear equations with a small number of unknowns.

If the time signals are truncated during data acquisition when the signal values have not reached a very small value, the unmeasured “tails” must be extrapolated and appended to the measured signals before the Fourier transform is taken. However, data extrapolation and interpolation may introduce errors. For irreversible tracers, the activity curve does not decay to zero at all while the input function decays to zero. Our proposed method does not apply and should be modified by using the derivatives of the time-activity-curves to replace the time-activity-curves.

To compare the standard non-linear iterative curve-fitting technique and the closed-form linear estimation methods, the iterative curve-fitting technique is more robust in a noisy environment but is strongly dependent on the initial conditions. Large bias on the estimated kinetic parameters (i.e., the k values) is common if the initial condition is not close enough to the true solution.

Many linear estimation methods are available. The linear coefficients are formed by the measured data and are thus influenced by noise. The linear estimation of the parameters is, in general, biased and noisy, even though no initial condition is needed.

For the two compartment model, the proposed method is sensitive to noise. Thus the result is not very accurate. The advantage of the proposed method is that the user does not need to specify the initial condition. On the other hand, the traditional curve-fitting method strongly depends on the initial condition. Therefore, the proposed method can be used to provide the initial condition that is close to the true solution, and then the traditional curve-fitting method is used to fine-tune the parameter estimation.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All authors read and approved the final manuscript.