Abstract
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 multicompartment 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 closedform formulas for estimating kinetic parameters for the twocompartment model. The proposed Fourier domain estimation method provides a unique solution, and offers very different noise response as compared to traditional nonlinear chisquared 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 Fourierdomain estimation method has closedform formulas. The proposed Fourierdomain curvefitting method does not require an initial condition, it minimizes a quadratic objective function and a closedform solution can be obtained. The noise is easier to control, simply by discarding the high frequency components, and emphasizing the DC component.
Keywords:
Kinetic parameter estimation; Dynamic imaging; Leastsquares estimation; Nuclear medicine imaging; Compartment modeling; Fourier transformIntroduction
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 [14], 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 multicompartment models is classic nonlinear leastsquares estimation [5]. Here, a leastsquares or weighted leastsquares objective function is iteratively minimized to obtain bestfit kinetic parameter estimates. Numerous curvefitting algorithms have been investigated for this application, including derivativebased or downhill simplex methods [5], ridgeregression [611], 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 leastsquares objective function for multicompartment models with noisy data is illformed, containing broad shallow valleys and (potentially) local minima. As such, careful implementation of the curvefitting algorithm with extensive iteration and handling of local minima is required to confidently find the global minimum. This results in computationaldemands that, while reasonable for fitting individual timeactivity curves, may become impractical for voxelwise parametric imaging where millions of fits need to be performed.
The conventional leastsquares curvefitting method is to match the measured timeactivitycurve with the calculated solution of the ordinary differential equations. The solution is in the form of a nonlinear function of time and the kinetic parameters. It is neither in the ODE form or the integral form. Much of the difficulty in the curvefitting 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 [1416], linear leastsquares (LLS) and generalized linear leastsquares (GLLS) methods [1720]. 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 timeactivity curves and blood input functions. Though the image noise is not correlated between timeframes, integrating the timeactivity curves in time gives equation errors that are not statistically independent, and that can bias the results (e.g. LLS); this correlationinduced bias can be reduced or eliminated by the iterative autoregressive filtering technique of GLLS. More recently, approaches employing temporal basis functions to reduce noise and improve parameter estimation have also been studied [2126]. Performance comparisons of these methods have been performed by Feng et al.[27] and more recently by Dai et al.[28]. Overall, nonlinear leastsquares 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 closedform formulas for estimating kinetic parameters for the twotissue 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 closedform solutions in the time domain [30], where the continuoustime differential equation is transformed into a discretetime difference equation by using the exact continuoustime expression of the compartment activity. This timedomain approach is fundamentally different from timedomain curvefitting approaches, and thus will likely have very different characteristics in terms of accuracy, speed, noise propagation and precision.
This paper introduces Fourierdomain closedform solutions for the estimation of kinetic parameters for the twotissue compartment model with 4 rate parameters. The similarity of the timedomain method and this Fourierdomain method are in the closedform 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 highfrequency noisy components can be readily discarded, and exponential functions never appear.
Methods and results
This section presents the standard method of nonlinear fitting to estimate the kinetic parameters, points out its potential deficiencies of the dependency on the initial conditions, and develops a closedform 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 curvefitting method
Figure 1 shows a generic twotissue compartment model with four rate constants: K_{1}, k_{2}, k_{3} and k_{4}. The model can be described by two firstorder differential equations:
Figure 1. A general twocompartmentmodel.
The above two equations are referred to as the state equations in system theory. The timeactivity curve for the blood input is B(t). There are two tissue compartments C_{1}(t) and C_{2}(t); these two compartments are not individually accessible and their timeactivity 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 timeactivitycurve C_{measured}(t) and the estimated timeactivitycurve C_{model}(t) evaluated by (6).
More precisely, , where C_{measured}(t_{i}) is the measured activity count at time t_{i} and w_{i} is the associated reciprocal variance. Also note C_{model}(t_{i}) need not be in the form of a convolution integral. We can alternately use an ODEsolver in conjunction with curvefitting to compute C_{model}(t_{i}) directly from the (1), (2), and (3) as needed in each curvefitting iteration.
Since the gradient of the objective function F in (11) is nonlinear with respect to the unknown parameters: K_{1}, k_{2}, k_{3} and k_{4}, 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 curvefitting method
In the following computer simulations, the blood input functions were in the form of [31]:
where A_{1} = 22.24, A_{2} = 8.36, A_{3} = 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 nonuniformly 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 60second 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 timeactivitycurve, 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 timeactivitycurve. 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 listmode, in which data are acquired continuously with a nondiscretized time stamp for each event. After the listmode 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 K_{1} = 0.4 min^{1}, k_{2} = 0.3 min^{1}, k_{3} = 0.2 min^{1} and k_{4} = 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, K_{1} can be anywhere from 0.1 to1.2, and k_{2} can be from 0.1 to 0.8, and so on. The values of k_{3} and k_{4} 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 T_{1/2} (=110 min) was the halflife 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 preselected nonuniformly distributed time instances, and the nonuniformly distributed blood and tissue timeactivitycurves are obtained.
Figure 2. Timeintegrated timeactivity curves C(t) 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 leastsquares fitting method was applied to the model (6) using MATLAB® and used the builtin function “nlinfit.” The MATLAB’s builtin function nlinfit allows negative outcomes. MATLAB also has another builtin function lsqnonlin that enforces the nonnegativity constraint on the outcomes and is widely used for nonnegative 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 closedform 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 curvefitting method with a good and a bad initial condition (using 1 noise realization for each case)
Table 2. The results of a standard iterative curvefitting method with a good and a bad initial condition (using 250 noise realizations for each case)
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 shorthand notations:
Using these shorthand notations, Eq. (14) becomes
Substituting (15), namely,, into (18) eliminates x_{4}, 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 innerproduct is defined by
A closedform solution can readily be obtained as
Finally, the twocompartment model kinetic parameters are estimated as
TAC extrapolation
If a prolonged data acquisition is not practical, the unmeasured “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_{end} is the time when data acquisition stops and λ_{n} (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_{end}) was the average value of the last 3 measurements.
Implementation
The computer procedure is discussed in this section. For the twocompartment model, the steps are:
Step 1: Acquire the blood input function B(t_{n}) and the compartment timeactivity curve C(t_{n}), for n = 0, 1, 2, …, N_{sample}1. Here the sampling interval can be nonuniform and N_{sample} is the total number of samples.
Step 2: Extrapolate the timeactivitycurve 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 4: Use a DFT (Discrete Fourier Transform) or FFT (Fast Fourier Transform) computer routine to calculate the DFT functions 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_{c} , corresponding to the continuous cutoff frequency ω_{c}.
Step 7: Form the matrices A and P as shown in (22) and (24).
Step 8: Solve for X using (26).
Step 9: Finally, use (27) to obtain estimated kinetic parameters K_{1} ~ k_{4}.
The proposed Fourier domain estimation method is in closedform 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 curvefitting 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_{c} was chosen to be 10. The nonuniformly sampled, onehour timeactivitycurves were extrapolated into twohour curves, and then linearly interpolated with a constant resampling interval of 5 seconds.
Some macroparameters have special clinical significance. For example, VD indicates volumeofdistribution, (K_{1}k_{3})/(k_{2} + k_{3}) reflects the net uptake, and (K_{1}k_{3})/( k_{2}k_{4}) or k_{3}/ k_{4} are related to the binding potential. Therefore, the estimation of some macroparameters 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 Fourierdomain closedform method. The requirement is that the Fourier integration interval must contain all nonzero activities. Therefore, it is expected that the traditional iterative nonlinear 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 leastsquares fitting usually gives satisfactory results. However, if an inappropriate initial condition is provided, the algorithm may converge to a wrong solution. When the closedform 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 closedform formula are discussed below.
Data integration effects
Averaging of the timeactivity 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 boxcar window function or unitstep 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 cutoff frequency
The selection of the cutoff s based on trialanderror, and is dependent on noise level. Ideally speaking, the cutoff 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 cutoff frequency is selected too high, the estimation will be affected by the noise fluctuation, resulting in different solutions with different noise realizations. If the cutoff 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 cutoff frequency selection is illustrated by the examples in Table 5, where a twocompartment model is considered with noise level α = 0.4.
Table 5. Estimation results for the twocompartment model parameters K_{1}, k_{2}, k_{3}, and k_{4}with noise level α = 0.4, using different cutoff frequencies and 250 noise realizations
Extension to onecompartment or multicompartment models
In general, if we have N, where N can be 1, compartments for the tissue model, we have a system of N firstorder 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 lefthandside 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 x_{1} … x_{N+N} 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 N1. The unknown coefficients are obtained by solving a leastsquares problem to approximate (35). Finally the kinetic parameters are estimated using relations between the coefficients x_{1} … x_{N+N} and the kinetic parameters. The proposed method becomes less effective as N increases, partially due to the nonlinear relationship between the k parameters and x parameters. For a large N, there may not be a closedform 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 nonidentifiable 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_{v} is the fraction of the contamination. Using the twocompartment 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_{v} as another unknown parameter, closedform formulas for the kinetic parameters as well as f_{v} can be obtained.
Relation to other linear methods
Generally speaking, a frequencydomain method is equivalent to a timedomain method. Using a cutoff in the frequencydomain is equivalent to performing smoothing in the timedomain. 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 zeromean, the DCcomponent is almost unaffected; 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 noiseheavilyaffected components (i.e., the high frequency components). We then weigh the remaining frequency components evenly, as an unweighted leastsquare objective function. We control the noise by dividing the frequency components in three categories: Trustworthy (i.e., the DC), somewhat trustworthy (i.e., the lowfrequencies), and nottrustworthy (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 closedform 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 lowfrequency 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
Timedomain curvefitting is the current stateoftheart in nuclear medicine kinetic estimation. Due to the nonlinear exponential functions, this curvefitting is sensitive to initial conditions (i.e., the initial solutions) for multicompartment model parameter estimation problems. The initial condition may make the algorithm converge to a wrong solution. In this paper, a Fourierdomain kinetic estimation method is proposed. The proposed Fourierdomain curvefitting method does not require an initial condition, it minimizes a quadratic objective function and a closedform solution can be obtained. The noise is easier to control, simply by discarding the high frequency components, and emphasizing the DC component. The proposed Fourierdomain estimation method has closedform 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 timeactivitycurves to replace the timeactivitycurves.
To compare the standard nonlinear iterative curvefitting technique and the closedform linear estimation methods, the iterative curvefitting 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 curvefitting 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 curvefitting method is used to finetune the parameter estimation.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors read and approved the final manuscript.
Acknowledgement
This work was supported in part by the Ben B. and Iris M. Margolis Foundation, NIH grants R01 HL108350, R01 CA135556, R01 HL50663, R01 EB007219 and by the Director, Office of Science, Office of Biological and Environmental Research of the US Department of Energy under contract DEAC0205CH11231. We thank Dr. Roy Rowley of the University of Utah for English editing.
References

Cherry SR, Sorenson JA, Phelps ME: Physics in Nuclear Medicine. 3rd edition. Philadelphia: Saunders; 2003.

Phelps ME: PET: Molecular Imaging and Its Biological Applications. New York: Springer Science; 2004.

Watabe H, Ikoma Y, Kimura Y, Naganawa M, Shidahara M: PET kinetic analysis—compartmental model.
Ann Nucl Med 2006, 20:583588. PubMed Abstract  Publisher Full Text

Gullberg GT, Reutter BW, Sitek A, Maltz J, Budinger TF: Dynamic single photon emission computed tomography — basic principles and cardiac applications.
Phys Med Biol 2010, 55:R111R191. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Press WH, Flannery BP, Teukolsky SA, Vetterling WT: Numerical Recipes in C. Cambridge: Cambridge University Press; 1988.

Byrtek M, O'Sullivan F, Muzi M, Spence AM: An adaptation of ridge regression for improved estimation of kinetic model parameters from PET studies.
IEEE 2003 Nuclear Science Symposium Conference Record 2003, 31203124.

Byrtek M, O'Sullivan F, Muzi M, Spence AM: An adaptation of ridge regression for improved estimation of kinetic model parameters from PET studies.

O'Sullivan F, Saha A: Use of ridge regression for improved estimation of kinetic constants from PET data.
IEEE Trans Med Imaging 1999, 18:115125. PubMed Abstract  Publisher Full Text

Yun Z, SungCheng H, Bergsneider M: Linear ridge regression with spatial constraint for generation of parametric images in dynamic positron emission tomography studies.
IEEE Trans Nucl Sci 2001, 48:125130. Publisher Full Text

Zhou Y, Huang SC, Bergsneider M, Wong DF: Improved parametric image generation using spatialtemporal analysis of dynamic PET studies.
Neuroimage 2002, 15:697707. PubMed Abstract  Publisher Full Text

Zhou Y, Endres CJ, Brasic JR, Huang SC, Wong DF: Linear regression with spatial constraint to generate parametric images of ligandreceptor dynamic PET studies with a simplified reference tissue model.
Neuroimage 2003, 18:975989. PubMed Abstract  Publisher Full Text

Wong KP, Meikle SR, Feng D, Fulham MJ: Estimation of input function and kinetic parameters using simulated annealing: application in a flow model.
IEEE Trans Nucl Sci 2002, 49:707713. Publisher Full Text

Yaqub M, Boellaard R, Kropholler MA, Lammertsma AA: Optimization algorithms and weighting factors for analysis of dynamic PET studies.
Phys Med Biol 2006, 51:42174232. PubMed Abstract  Publisher Full Text

Blomqvist G: On the construction of functional maps in positron emission tomography.
J Cereb Blood Flow Metab 1984, 4:629632. PubMed Abstract  Publisher Full Text

Carson RE, Huang SC, Green ME: Weighted integration method for local cerebral blood flow measurements with positron emission tomography.
J Cereb Blood Flow Metab 1986, 6:245258. PubMed Abstract  Publisher Full Text

Yokoi T, Kanno I, Iida H, Miura S, Uemura K: A new approach of weighted integration technique based on accumulated images using dynamic PET and H2(15)O.
J Cereb Blood Flow Metab 1991, 11:492501. PubMed Abstract  Publisher Full Text

Chen K, Lawson M, Reiman E, Cooper A, Feng D, Huang SC, Bandy D, Ho D, Yun LS, Palant A: Generalized linear least squares method for fast generation of myocardial blood flow parametric images with N13 ammonia PET.
IEEE Trans Med Imaging 1998, 17:236243. PubMed Abstract  Publisher Full Text

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:3143. PubMed Abstract  Publisher Full Text

Ho D, Feng D: Rapid algorithms for the construction of cerebral blood flow and oxygen utilization images with oxygen15 and dynamic positron emission tomography.
Comput Methods Programs Biomed 1999, 58:99117. PubMed Abstract  Publisher Full Text

Wen L, Eberl S, Fulham MJ, Feng DD, Bai J: Constructing reliable parametric images using enhanced GLLS for dynamic SPECT.
IEEE Trans Biomed Eng 2009, 56:11171126. PubMed Abstract  Publisher Full Text

Boellaard R, Knaapen P, Rijbroek A, Luurtsema GJ, Lammertsma AA: Evaluation of basis function and linear least squares methods for generating parametric blood flow images using 15Owater and Positron Emission Tomography.
Mol Imaging Biol 2005, 7:273285. PubMed Abstract  Publisher Full Text

Gunn RN, Gunn SR, Turkheimer FE, Aston JA, Cunningham VJ: Positron emission tomography compartmental models: a basis pursuit strategy for kinetic modeling.
J Cereb Blood Flow Metab 2002, 22:14251439. PubMed Abstract  Publisher Full Text

Hong Y: T and Fryer T D: Kinetic modelling using basis functions derived from twotissue compartmental models with a plasma input function: general principle and application to [18 F]fluorodeoxyglucose positron emission tomography.
Neuroimage 2010, 51:164172. PubMed Abstract  Publisher Full Text

Reader AJ, Sureau FC, Comtat C, Trebossen R, Buvat I: Joint estimation of dynamic PET images and temporal basis functions using fully 4D MLEM.
Phys Med Biol 2006, 51:54555474. PubMed Abstract  Publisher Full Text

Verhaeghe J, Van de Ville D, Khalidov I, D'Asseler Y, Lemahieu I, Unser M: Dynamic PET reconstruction using wavelet regularization with adapted basis functions.
IEEE Trans Med Imaging 2008, 27:943959. PubMed Abstract  Publisher Full Text

Watabe H, Jino H, Kawachi N, Teramoto N, Hayashi T, Ohta Y, Iida H: Parametric imaging of myocardial blood flow with 15Owater and PET using the basis function method.
J Nucl Med 2005, 46:12191224. PubMed Abstract  Publisher Full Text

Feng D, Ho D, Chen K, Wu LC, Wang JK, Liu RS, Yeh SH: An evaluation of the algorithms for determining local cerebral metabolic rates of glucose using positron emission tomography dynamic data.
IEEE Trans Med Imaging 1995, 14:697710. PubMed Abstract  Publisher Full Text

Dai X, Chen Z, Tian J: Performance evaluation of kinetic parameter estimation methods in dynamic FDGPET studies.
Nucl Med Commun 2011, 32:416. PubMed Abstract  Publisher Full Text

Zeng GL, Kadrmas DJ, Gullberg GT: Fourier domain closedform formulas for estimation of kinetic parameters in multicompartment models.
2011 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC 2011)32093216.

Zeng GL, Gullberg GT, Kadrmas DJ: Closedform kinetic parameter estimation solution to the truncated data problem.
Phys Med Biol 2010, 55:74537468. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Feng D, Huang SC, Wang X: Models for computer simulation studies of input functions for tracer kinetic modeling with positron emission tomography.
Int J Biomed Comput 1993, 32:95110. PubMed Abstract  Publisher Full Text

Oriuchi N, Tomiyoshi K, Ahmed K, Sarwar M, Tokunaga M, Suzuki H, Watanabe N, Hirano T, Shibasaki T, Tamura M, Endo K: Independent Thallium201 accumulation and Fluorine18Fluorodeoxyglucose metabolism in glioma.
J Nucl Med 1996, 37:457462. PubMed Abstract  Publisher Full Text

Ichise M, Toyama H, Innis RB, Carson RE: Strategies to improve neureceptor parameter estimation by linear regression analysis.
J Cereb Blood Flow Metab 2002, 22:12711281. PubMed Abstract  Publisher Full Text

Gunn RN, Gunn SR, Cunningham VJ: Positron emission tomography compartmental models.
J Cereb Blood Flow Metab 2001, 21:635652. PubMed Abstract  Publisher Full Text

Franklin GF, Powell JD, EmamiNaeini A: Feedback Control of Dynamic Systems. 4th edition. Upper Saddle River, New Jersey: PrenticeHall Inc; 2002.