Estimation of input functions from dynamic [18F]FLT PET studies of the head and neck with correction for partial volume effects

Background We present a method for extracting arterial input functions from dynamic [18F]FLT PET images of the head and neck, directly accounting for the partial volume effect. The method uses two blood samples, for which the optimum collection times are assessed. Methods Six datasets comprising dynamic PET images, co-registered computed tomography (CT) scans and blood-sampled input functions were collected from four patients with head and neck tumours. In each PET image set, a region was identified that comprised the carotid artery (outlined on CT images) and surrounding tissue within the voxels containing the artery. The time course of activity in the region was modelled as the sum of the blood-sampled input function and a compartmental model of tracer uptake in the surrounding tissue. The time course of arterial activity was described by a mathematical function with seven parameters. The parameters of the function and the compartmental model were simultaneously estimated, aiming to achieve the best match between the modelled and imaged time course of regional activity and the best match of the estimated blood activity to between 0 and 3 samples. The normalised root-mean-square (RMSnorm) differences and errors in areas under the curves (AUCs) between the measured and estimated input functions were assessed. Results A one-compartment model of tracer movement to and from the artery best described uptake in the tissue surrounding the artery, so the final model of the input function and tissue kinetics has nine parameters to be estimated. The estimated and blood-sampled input functions agreed well when two blood samples, obtained at times between 2 and 8 min and between 8 and 60 min, were used in the estimation process (RMSnorm values of 1.1 ± 0.5 and AUC errors for the peak and tail region of the curves of 15% ± 9% and 10% ± 8%, respectively). A third blood sample did not significantly improve the accuracy of the estimated input functions. Conclusions Input functions for FLT-PET studies of the head and neck can be estimated well using a one-compartment model of tracer movement and TWO blood samples obtained after the peak in arterial activity.


Background
The role of positron emission tomography (PET) in cancer staging and assessment of treatment response continues to grow [1]. The radiotracer 3 -deoxy-3 -[ 18 F]fluorothymidine (FLT) is a marker for cell proliferation as it is phosphorylated by thymidine kinase-1, which is selectively expressed in the S, G2 and M phases of the cell http://www.ejnmmires.com/content/3/1/84 to patients [5]. Consequently, other sources of arterial input functions have been explored, such as TACs taken from the vicinity of large blood vessels in dynamic PET images [6,7]. However, accurate determination of input functions from PET images is hindered by the 'partial volume effect' (PVE), a term that refers to two distinct but related phenomena: spillover and the tissue-fraction effect [8]. Spillover is caused by the finite spatial resolution of the imaging system and manifests as a blurring of the imaged distribution of activity. The tissue-fraction effect arises because the activity reported in each PET voxel is an average of the corresponding volume, which usually comprises tissues of different types and activity concentrations.
The diameter of the largest blood vessels in the head and neck is around 5 mm [9], which is comparable to the resolution of PET scanners and the dimensions of PET voxels. Some methods used to extract input functions from PET images attempt to explicitly correct the PVE using anatomical information taken from computed tomography (CT) or magnetic resonance scans [10][11][12]. Other approaches, most commonly used in brain imaging, assume that the activity measured in a region represents a mixture of blood and surrounding tissue activities and implicitly correct the PVE by modelling the relative contributions of blood and tissue to the total imaged regional activity [13][14][15][16].
Models of both the input function and tracer movement between the blood vessel and the surrounding tissue can be simultaneously fitted to a total regional TAC taken from the vicinity of the blood vessel in the dynamic PET images. This technique has shown promise for the radiotracer 2-deoxy-2-[ 18 F]fluoro-D-glucose (FDG) in human [13] and animal [15] studies, and for the radiotracers FLT [16] and 18 F-fluoromisonidazole (Fmiso) [17] in human studies. In this approach, the input function is represented as a mathematical formula, and activity in the surrounding tissue is linked to it via a compartment model. The parameters of the input formula and compartment model are estimated simultaneously so that the combined modelled activity best matches the regional TAC obtained from the dynamic PET images. A challenge for simultaneous estimation is that the number of parameters may be large and a unique solution not identifiable. Further, the previous studies have used a variety of compartmental (and empirical) models to describe tracer movement between a vessel and the surrounding tissue. An optimal model, which best balances accurate description of the tissue TAC and parsimony of parameters, has not yet been established.
This work identifies an appropriate model of tracer movement, then uses the model to simultaneously estimate arterial input functions for bolus injections of FLT from TACs of carotid artery regions drawn on dynamic FLT-PET images. We also explore the use of small numbers of arterial concentration measurements obtained from blood samples to refine possible solutions and thus improve the precision of the resulting arterial TAC estimates. Optimal times for blood sampling have been assessed.

Methods
Written informed consent was obtained from all participants, and the study was approved by the local Research Ethics Committee and UK Administrat0ion of Radioactive Substances Advisory Committee.

Patient details and image acquisition
Six datasets were collected for four patients diagnosed with stage II-III head and neck squamous cell carcinoma. Patients were scanned at the PET imaging centre of St Thomas' Hospital, London using a GE Healthcare Discovery VCT PET-CT scanner (Waukesha, WI, USA). After CT scanning of the head and neck region, patients were injected with 2.59 MBq/kg of [ 18 F]FLT (up to a maximum of 185 MBq), and PET imaging begun. List-mode PET data were binned into 24 frames (8 × 15 s, 4 × 30 s, 6 × 60 s, 2 × 300 s and 4 × 600 s) for four of the dynamic image sets, while for two image sets, data were binned into 103 frames (30 × 2 s, 12 × 10 s, 6 × 20 s and 55 × 60 s).
Images were reconstructed using Fourier rebinning, and a two-dimensional filtered back projection algorithm to ensure quantitative accuracy [18]. Corrections for attenuation, scatter and dead time were performed during sinogram histogramming and reconstruction.
An image array of 128 × 128 × 47 voxels per frame was used, with dimensions of 5.47 × 5.47 × 3.27 mm 3 . The transaxial full-width half-maximum (FWHM) and full-width tenth-maximum (FWTM) values of the scanner at 1 cm along the transaxial axis were 4.9 and 9.5 mm, respectively, and the axial FWHM and FWTM were 4.8 and 10.6 mm, respectively. The voxel spacing of the coregistered CT images was 0.98 × 0.98 × 3.27 mm 3 .
Blood samples were collected through a cannula placed in the radial artery, and 18 F activity was measured using an Allogg Arterial Blood Sampling System (Allogg AB, Mariefred, Sweden). For the first 24 min, samples were taken at a rate of 1/s. A further two samples were taken at 38 and 56 min post-injection. Blood-sampled input functions, b(t), were constructed from 18 F concentrations in the samples.

Modelling the imaged TAC of the arterial region using the blood-sampled input function
For each dataset, the carotid artery was contoured on eight contiguous CT slices and the contours are transferred to the PET images. An example of the outlined carotid artery is shown in Figure 1. The corresponding PET voxels, also shown in Figure 1, that included the http://www.ejnmmires.com/content/3/1/84 artery also comprised a significant fractional volume of extravascular tissue. The arterial voxels were identified on each of the eight slices, and the fractional volumes of artery and surrounding tissue, v a and v tis , respectively, were determined. The carotid artery typically constituted 30% to 50% of the arterial voxel volume in the images used in this study. The total activity of the arterial voxels, m tot , can be modelled as where C a and C tis are models of the input function and the TAC of the surrounding tissue, respectively, and f is a compartmental (or empirical) model of the tissue TAC with input function C a . However, the total regional activity recorded in the delineated PET voxels also comprises spill-in from the tissue in the next layer of voxels out, less the spill-out from the arterial voxels. Thus, the total activity also needs to be modified to correct for spillover, as described in the 'Modelling of spillover' section: We have tested four different compartmental models: (a) irreversible uptake of tracer in a single non-vascular compartment, (b) reversible tissue uptake, (c) reversible uptake of tracer in the first non-vascular compartment from which tracer was irreversibly taken up in a second non-vascular compartment and (d) reversible uptake of tracer in the first and second non-vascular compartments. These models have one, two, four and five parameters, respectively. We have also tested a two-parameter, phenomenological model of tracer movement used by Backes et al. [16]. For each dataset, the five m tot models were fitted to the mean activity, PET tot , recorded in the arterial voxels using the corresponding blood-sampled data as the input function for the compartmental model. The MATLAB function 'lsqnonlin' (The MathWorks, Natick, MA, USA) was used to fit m tot to PET tot by minimizing the objective function where t i is the time at the midpoint of the ith of N measurement intervals, and each weight w i is inversely proportional to the variance of the imaged tracer concentration averaged over the ith measurement interval.
The optimal model was identified as that with the lowest Akaike information criterion corrected for small sample size (AIC c ) [19]. The AIC c values for all datasets are shown in Table 1. For all datasets, model B had the lowest AIC c value, indicating that a compartmental model that accounted for tracer movement to and from the vessel was necessary and sufficiently complex to describe the data. The equation of state for model B is where K 1 and k 2 are kinetics parameters describing the rate of tracer movement to and from the tissue, respectively.

Simultaneous estimation of the input function and tissue TAC from PET tot
Having found the model m tot which provides the best estimate of the activity recorded in the arterial voxels, we again fitted this model to PET tot in order to extract the input function by minimizing the objective function given in Equation 2. In this process, it was assumed that the unknown input function is described by a mathematical model whose parameters are simultaneously estimated along with the kinetics parameters of the compartmental model. The input function C a can conveniently http://www.ejnmmires.com/content/3/1/84 and accurately be represented by the seven-parameter function [20]: An example of a blood-sampled input curve is shown in Figure 2, together with a fit of Equation 4. Simpler mono-and bi-exponential models of the input function (both with a time-delay parameter for tracer delivery) were also tested but could not adequately describe the blood-sampled input functions. Thus, the final model of the total activity of the arterial pixels has nine parameters to be estimated: two kinetic parameters and seven input function parameters.
Blood samples obtained at a single time point were incorporated into the process of fitting m tot to PET tot by where b (t k ) is the measured arterial concentration at time t k andw is the mean value of the weighting factors w i . The constant α is a ratio by which the squared difference between the estimated and measured blood concentrations has more influence on the fitting process than the average of the squared residual errors of the modelled m tot (t) values. A range of α values from 2 to 11 was considered and changes in the accuracy of the C a models as a function of α were assessed using analysis of variance (ANOVA). The nine independent parameters of the full m tot model are reduced to eight if the model is constrained to match a blood sample value at a particular time point. Thus, in order to avoid under-or over-emphasizing the fit at that time point, a reasonable α value would be around N × 1/8 (that is, 3 or 13 for N = 24 or 103), so that the weight given to the fit of the blood-sample point is one eighth that of the total weight of the fit of the model to PET tot (t). For additional blood samples, fitting was achieved by adding further terms to Equation 5, so that for three samples taken at times t k , t l and t m , the objective function was used. For two and three blood samples, reasonable α values would be 2N/ (2 × 7) and 3N/ (3 × 6) (3 and 4 respectively for N = 24, and 14 and 17 respectively for N = 103).

Initial values of model parameters used for simultaneous estimation
The initial values of the kinetics parameters K 1 and k 2 used in the simultaneous estimation were the median values of the parameters obtained when Equation 1 (using model B) was directly fitted to PET tot using b(t) as the http://www.ejnmmires.com/content/3/1/84 input function. The initial values used for the input function parameters were the medians of those obtained by fitting Equation 4 to the PET tot (t) data and are termed Init PET . We also studied the use of the median values Init BSIF , obtained from directly fitting Equation 4 to the blood-sampled input function, as initial values for the estimated input function.

Metrics used for evaluation of results
To evaluate agreement between a blood-sampled input function b(t) and a simultaneously estimated input function C a (t) we used a scale-independent measure of the root-mean-square (RMS) error Equation 4 to the bloodsampled input function b(t). The RMS norm parameter therefore represents the goodness-of-fit of the estimated input function to the blood-sampled data relative to the goodness-of-fit of a direct fit of Equation 4 to the bloodsampled data. Errors on areas under the curve (AUC) of the estimated input functions were also calculated. The input functions were divided into an early region, characterized by a sharp peak in activity, and a tail during which tracer concentrations in the blood and tissue gradually equilibrated. The boundary between the two regions was defined by the maximum value of the second derivative of b(t), typically around 6 min. The AUC error for each region was expressed as the modulus of the difference between the areas under the estimated and blood-sampled input functions, calculated as a percentage of the area under the blood-sampled input function.

Modelling of spillover
Spillover was incorporated into the models of the total regional activity using the geometric transfer matrix (GTM) method [10]. This approach assumes that the activity is distributed over a number of regions, each with homogeneous concentration of activity. The GTM relates a vector of the true concentrations in the regions at time point t k , C true (t k ), to a vector of measured concentrations C meas (t k ), via where the elements g ij of the matrix G represent the contribution of the true activity in region i to the activity recorded in region j. The GTM is usually inverted to recover the true concentrations in the regions of interest: The activity recorded in the layer of voxels surrounding the arterial voxels, termed the external voxels, was significantly lower than the activity of the non-vascular tissue within the arterial voxels, termed the surrounding tissue, particularly around the peak in arterial activity concentration, as illustrated by Additional file 1: Figure S1. Thus, three regions of interest were considered for the GTM method: the artery, the surrounding tissue and the external voxels. Anatomical information from the CT images was used to construct a high-resolution geometric model of each region of interest. The activity in each specific region was sequentially set to 1, and all other regions were set to 0. For each region i, this model was then convolved with a model of the point-spread function (PSF) of the PET scanner, and the intersection fraction of region i with each region j was calculated to generate the elements g ij .
The simultaneous estimation process aims to estimate the input function by fitting the model of Equation 1 to the total activity recorded in the delineated arterial voxels, which is affected by spillover. Equation 6 was therefore used to simulate the effects of spillover on the 'true' estimated concentrations of the artery and surrounding tissue concentrations before matching the total activity of these regions, m tot (t k ), to that of the imaged region, PET tot (t k ).
As the activity recorded in the arterial voxels also includes spill-out from the external voxels, the true concentration of activity in these voxels, C ext (t k ), is required to calculate m tot (t k ). C ext (t k ) was estimated using Equation 7 to correct the activity measured in these voxels for spill-out from the arterial voxels. For this purpose, the activity in both the artery and surrounding tissue was set as the mean of the total activity in the arterial voxels. A vector of the estimated values of the arterial activity, surrounding tissue activity and the true activity in the external voxels was then multiplied by G to obtain a model of m tot (t k ) with simulation of spillover. Models of m tot (t k ) without simulation of spillover were also generated to examine the importance of explicitly simulating the effects of spillover on the accuracy of the estimated input functions.
The effects of the reconstruction algorithm and associated filters were modelled as a composite of two Gaussian functions fitted to the PSF of the PET scanner. This kernel was assumed to be stationary as all regions of interest lay within 5 cm of the centre of the transaxial field of view.

Statistical analysis
ANOVA was used to study the influence of α on RMS norm values of simultaneously estimated input functions. For simultaneous estimation using a single blood sample, input functions were generated for different sampling times t k , and median values and interquartile ranges of their RMS norm values and AUC errors were tabulated, http://www.ejnmmires.com/content/3/1/84  splitting the sampling times into 'early' (0< t ≤ 2 min), 'intermediate' (2< t ≤ 8 min) and 'late' (8< t ≤60 min) intervals. The interquartile range describes the change in a metric with sampling time t k throughout an interval and also reflects the variability introduced by noise on the measured arterial tracer concentration.
The F test was used to compare residual errors for input functions generated using one versus two and two versus three blood samples. The appropriate (combination of ) sampling times were identified as those times that resulted in the lowest median RMS norm values for all six datasets. When two blood samples were used, median and interquartile ranges of RMS norm and AUC errors were tabulated for the six possible combinations of early, intermediate and late blood-sampling intervals, and ANOVA was used to assess differences in these metrics between different combinations of sampling intervals. RMS norm distributions of input functions generated with and without simulating spillover were also compared using ANOVA. Table 2 summarises the metrics of agreement between the blood-sampled and estimated input functions when blood samples were taken at (combinations of ) times found to minimise RMS norm values. Good agreement was achieved for all datasets (RMS norm values of 1.0 to 1.5, median value 1.1) when arterial concentrations measured at two time points were used in the estimation process, but inclusion of a third sample did not significantly improve agreement for five of the six datasets (one p value of 0.002, all other p values > 0.35). However, the input functions estimated without using any arterial concentration measurements showed notably worse agreement with the blood-sampled input functions (median RMS norm 2.9) and the use of one measurement was not sufficient to achieve reasonable agreement (median RMS norm 2.0). The F tests showed that inclusion of a second sample significantly improved the accuracy of five of the six estimated input functions datasets (p values < 0.007). Figure 3 shows an example of input functions simultaneously estimated for one of the datasets using two and three concentration measurements and default α values of N×1/8, 2/(2 × 7), 3/(3 × 6) respectively (initial results showed no significant dependence of RMS norm on α; all p values ≥ 0.35).

Simultaneous estimation of the arterial input function and surrounding tissue TAC
The accuracy of the input functions simultaneously estimated using zero or one arterial concentration measurements depends on the initial estimates of the model parameters; more accurate estimates were obtained by starting from Init BSIF than Init PET . When only one arterial concentration measurement was used in the simultaneous estimation process, more accurate estimates of the input function were generally obtained using a measurement made in the late time interval (8< t ≤ 60 min). However, the use of one concentration measurement often improved the accuracy of the estimated input function only in the vicinity of the time point for which the measurement was obtained, whilst the accuracy over other intervals was unchanged or even worsened. In contrast, the accuracy of input functions estimated using two or three arterial concentration measurements was independent of initial parameter values. The lowest RMS norm values were obtained using combinations of intermediate and late times of concentration measurement, as shown in Figure 4. Inclusion of the second concentration measurement taken in the intermediate interval improved the agreement between the bloodsampled and estimated input functions in the peak region without reducing the agreement in the tail region.

Spillover
For one, dataset, explicitly simulating spillover improved the agreement between input functions estimated with two blood samples ( p value = 0.03), but distributions of RMS norm values were not significantly different for input functions estimated with and without spillover for the other five datasets. Metrics for agreement between blood-sampled and simultaneously estimated input functions obtained without simulation of spillover are summarized in Table 3. The (combinations of ) sampling times are those used to generate the results presented in Table 2.

Discussion
It is clearly advantageous to generate input functions through non-invasive means, but the partial volume effect prevents direct extraction of accurate input functions from PET images. There is often a substantial tracer uptake in the tissue surrounding the arterial wall, and so any correction strategy for the PVE that ignores tracer uptake in the surrounding tissue is unlikely to accurately determine the input function. The method developed here enables accurate and robust estimation of input functions using image data and two blood samples. A one-compartment model that accounts for movement of the tracer to and from the artery is complex Table 3 Metrics for agreement of simultaneously estimated and blood-sampled input functions without simulation of spillover  (15) 18 (14) 15 Median (median interquartile range) of RMS norm values and errors of the area under the peak and tail of the estimated input function, as the moduli of the percentages of the areas under the blood-sampled input function (calculated without simulation of spillover). http://www.ejnmmires.com/content/3/1/84 enough to describe the tracer kinetics but does not include unnecessary parameters that might hinder the estimation process. The complete model of the input function and tissue kinetics has nine parameters, seven for the input function model and two for the kinetics parameters, to be simultaneously estimated from the TAC of voxels that include the artery and two tracer concentrations obtained from blood samples taken at different time points. When two arterial tracer concentration measurements, obtained at suitable times, are used in the simultaneous estimation process, the accuracy of the estimated input function is independent of the initial values chosen for the model parameters. If only one blood sample is used, the accuracy of the estimated input function depends on the choice of initial values, and the blood sample usually improves the accuracy of the estimated input function only around the time point at which the sample was taken. A combination of blood samples taken at intermediate and late times reduces the median RMS norm values of the estimated input functions to only 10% greater than those obtained by directly fitting the input function model (Equation 4) to the arterial TACs obtained by continual blood sampling. Explicitly modelling the effects of spillover does not generally improve the accuracy of the input functions estimated using the proposed method.
The model described here works well for FLT PET headand-neck images, from which TACs can be obtained for regions containing the carotid artery. However, it may not be appropriate for other tracers and/or imaged regions [12,21] for which the kinetics of tracer transport may differ. The process of model validation and selection, and exploration of the influence of direct measurements of blood concentrations on the accuracy of the estimated input functions, should be repeated for other radiotracers and imaging sites.