Journal of Cerebral Blood Flow & Metabolism (2015) 35, 1368–1379 © 2015 ISCBFM All rights reserved 0271-678X/15 $32.00 www.jcbfm.com
ORIGINAL ARTICLE
Model-free quantification of dynamic PET data using nonparametric deconvolution Francesca Zanderigo1, Ramin V Parsey2 and R Todd Ogden3,4,5 Dynamic positron emission tomography (PET) data are usually quantified using compartment models (CMs) or derived graphical approaches. Often, however, CMs either do not properly describe the tracer kinetics, or are not identifiable, leading to nonphysiologic estimates of the tracer binding. The PET data are modeled as the convolution of the metabolite-corrected input function and the tracer impulse response function (IRF) in the tissue. Using nonparametric deconvolution methods, it is possible to obtain model-free estimates of the IRF, from which functionals related to tracer volume of distribution and binding may be computed, but this approach has rarely been applied in PET. Here, we apply nonparametric deconvolution using singular value decomposition to simulated and test–retest clinical PET data with four reversible tracers well characterized by CMs ([11C]CUMI-101, [11C]DASB, [11C]PE2I, and [11C]WAY-100635), and systematically compare reproducibility, reliability, and identifiability of various IRF-derived functionals with that of traditional CMs outcomes. Results show that nonparametric deconvolution, completely free of any model assumptions, allows for estimates of tracer volume of distribution and binding that are very close to the estimates obtained with CMs and, in some cases, show better test–retest performance than CMs outcomes. Journal of Cerebral Blood Flow & Metabolism (2015) 35, 1368–1379; doi:10.1038/jcbfm.2015.65; published online 15 April 2015 Keywords: binding; full quantification; model-free; PET; singular value decomposition
INTRODUCTION Positron emission tomography (PET) is a nuclear imaging technology that uses radioactively labeled molecules (tracers) that are injected into the body and bind, although not exclusively, to a specific target to quantify in vivo biological processes, such as blood flow, metabolism, and the distribution of proteins in the brain and the body. PET is an invaluable tool to establish the genetic and cellular basis of diseases. Fully quantitative PET allows estimation of outcome measures such as, in the case of a reversible tracer, the total volume of distribution (VT) of the injected tracer and a series of binding potentials (BPs) between tracer and target, to indirectly determine the available amount of the target.1 These BPs are in various ways related to the density of available target (Bmax) through the affinity of the tracer for the available binding sites (1/KD).1 Quantification of these outcome measures relies on describing the PET signal associated with the tracer in the tissue (i.e., the time activity curve, TAC) as the convolution between the concentration of free (not protein bound) unmetabolized tracer in the arterial plasma (i.e., the metabolite-corrected input function) and a function specific to the tissue, the impulse response function (IRF).2 The metabolite-corrected input function represents how much of the injected tracer is actually available to bind with the target. IRF is estimated from measures of the input and the tissue PET signal, and would not otherwise be known. Traditionally, PET quantification methods have used parametric approaches, which require specifying a certain parametric form for the IRF (e.g., a mathematical function derived from a physiologic
model of the system, or that approximates the system behavior). This approach therefore relies on fairly strong assumptions, which are not always substantiated in reality, and a priori knowledge of the system under investigation. Traditional PET parametric analysis quantifies these outcome measures by interpreting the PET signal in the brain tissue using compartment models (CMs) to describe the system IRF (the kinetics) of the tracer, in particular a reversible or irreversible two-tissue compartment model (2TCM).1 The tracer signal in the brain tissue is modeled as the sum of two components: part of the tracer is specifically bound to the target of interest (specific binding); part is nonspecifically bound to macromolecules other than the target, or free in the tissue water (nonspecific binding).1,2 Identification of this model involves estimating four rate constants (K1, k2, k3, k4), which describe how quickly the tracer transits from ‘being in the plasma’, to ‘being free to bind in the tissue’, to ‘being specifically bound to the target’.1,2 However, these constants are not always all identifiable, and as a result directly estimating the specific binding component in the model can often be unstable. Therefore, full quantification of this model also requires identifying a reference region (RR) in the brain that is completely devoid of the target.1,2 In a proper (i.e., measurable, valid, reliable, independent of treatment effects and groups) RR devoid of the target of interest, the tracer can only be nonspecifically bound to macromolecules other than the target,1,2 and its kinetics can be described with a one-tissue compartment model1 (1TCM) regulated by only two rate constants (K1’ and k2’). The tracer volume of distribution in such a region is assumed to represent (and used to estimate) a common nonspecific binding
1 Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute and Columbia University, New York, New York, USA; 2Department of Psychiatry and Radiology, Stony Brook University, Stony Brook, New York, USA; 3Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute, New York, New York, USA; 4Department of Psychiatry, Columbia University, College of Physicians and Surgeons, New York, New York, USA and 5Department of Biostatistics, Columbia University, Mailman School of Public Health, New York, New York, USA. Correspondence: Dr F Zanderigo, Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute and Columbia University, 1051 Riverside Drive, New York, NY 10032, USA. E-mail: [email protected] or [email protected] Received 1 December 2014; revised 25 February 2015; accepted 18 March 2015; published online 15 April 2015
Model-free nonparametric PET quantification F Zanderigo et al
1369 across the brain, and helps isolate the component that is only specifically bound to the target. Often, however, CMs or graphical approaches derived from CMs2 either do not properly describe the tracer kinetics, or are not identifiable, leading to nonphysiologic estimates of the VT and of the binding to the target. On the other hand, nonparametric deconvolution approaches do not impose any specific formulation or model for the IRF, and only rely on very loose requirements, often matched in the case of physiologic systems, for the IRF (e.g., degree of smoothness; monotonicity). Nonparametric deconvolution approaches, therefore, allow obtaining model-free estimates of the IRF, from which functionals related to tracer VT and BPs can be computed. Nevertheless, they have been rarely adopted in PET. O'Sullivan et al3 proposed the use of nonparametric residue analysis and demonstrated it with simulations and data using the irreversible tracer [18F]FDG. Only recently, Jiang et al4 proposed a nonparametric approach based on functional principal component analysis and demonstrate it with simulations and data with the tracer [11C]-diprenorphine. Several approaches have been developed that fall between these two extremes.5–7 Generally, these involve expressing the IRF in terms of basis functions that, although derived from the principles of linear time-invariant systems and based on CMs,8 do not impose any predefined compartmental structure. Similar to what we propose using nonparametric deconvolution approaches, these basis function approaches allow for calculation of functionals that can be computed from the IRF as PET outcome measures. Here, we propose the use of nonparametric deconvolution via singular value decomposition (SVD)9,10 to compute functionals of the IRF of PET reversible tracers that are associated with tracer VT and BPs, and extensively compare accuracy, reproducibility, reliability, and identifiability of these different functionals with traditional outcome measures derived using CMs. SVD is an easy-to-implement nonparametric deconvolution approach that is widely used, for example, in the quantification of Dynamic Susceptibility Contrast-Magnetic Resonance Imaging data9,10 and has been used before to obtain regularized least-square solutions.7,11–13 We show relative advantages and disadvantages of using outcome measures derived from either nonparametric deconvolution or traditional parametric approaches, using simulated data and four test–retest clinical datasets with the tracers [11C]CUMI-101 (CUMI), [11C]DASB (DASB), [11C]PE2I (PE2I), and [11C] WAY-100635 (WAY). MATERIALS AND METHODS Subjects Four previously acquired test-retest datasets were considered: 12 scans with CUMI,14 20 with DASB,15 14 with PE2I,16 and 10 with WAY.17 The studies were performed in accordance with the Declaration of Helsinki and The Institutional Review Boards of Columbia University Medical Center and New York State Psychiatric Institute approved the protocol. Subjects gave written informed consent after receiving an explanation of the study.
Positron Emission Tomography Protocol Tracers preparation, emission data acquisition and reconstruction, and determination and fit of the metabolite-corrected input function were obtained as previously described.14–17
Image Analysis Images were analyzed using Matlab Release 2013b (The Mathworks, MA, USA) with extensions to FLIRT version 5.2 (ref. 18) and BET version 1.2 (ref. 19). Motion correction was applied. An RR and several anatomic target regions of interest (ROIs) were traced based on brain atlases and published reports:20,21 13 ROIs for CUMI (hippocampus, entorhinal cortex, insula, posterior parahippocampal gyrus, temporal lobe, amygdala, medial © 2015 ISCBFM
prefrontal cortex, cingulate, orbital, dorsolateral prefrontal cortex, parietal, raphe, occipital lobe; RR: cerebellar gray matter); 6 ROIs for DASB (midbrain, anterior cingulate, amygdala, dorsal putamen, hippocampus, thalamus; RR: cerebellar gray matter); 4 ROIs for PE2I (dorsal caudate, dorsal putamen, ventral striatum; RR: cerebellum); and 13 ROIs for WAY (raphe; anterior cingulate; amygdala; cingulate, dorsolateral prefrontal cortex, hippocampus, insula, medial prefrontal cortex, occipital lobe, orbital, parietal, posterior parahippocampal gyrus, temporal lobe; RR: cerebellar white matter). Selecting these ROIs restricted the analysis to regions for which the binding is expected to be relatively high (with the exception of the RR).
Simulated Data In addition to the clinical test–retest data sets, we considered simulations designed to generate data that imitate characteristics of real data in the case of 1TCM and 2TCM, respectively. For the 1TCM, we used the metabolite-corrected input function measured in one of the considered DASB scans, a tracer whose kinetics is well described by a 1TCM,15 and set the K1 and k2 values in six selected regions (cerebellar gray matter, temporal lobe, hippocampus, dorsal caudate, amygdala, and ventral striatum) based on existing data on DASB healthy control subjects.15 For the 2TCM, we used the metabolitecorrected input function measured in one of the considered CUMI scans, a tracer whose kinetics is well described by a 2TCM,14 and set the K1 to k4 values in five selected regions (cerebellar gray matter, hippocampus, temporal lobe, occipital lobe, and cingulate) based on existing data on CUMI healthy control subjects.14 In both cases, we simulated Gaussian noise with zero mean and standard deviation (SD) calculated as SD(ti) = α. σ/w(ti), where w(ti) is the root square of the duration of the time frame during which the PET data corresponding to time ti has been acquired, and σ indicates the level of noise that we derived from the weighted residual sum of square values obtained by fitting DASB and CUMI scans, respectively, with CMs. Increasing α allows us to study the various estimation techniques in situations with proportionately more noise than what we encounter in our real-data situations. We simulated five noise conditions setting α equal to 1, 2.5, 5, 7.5, and 10, respectively. In all cases and noise conditions, we simulated 1,000 TACs for each of the selected regions. The complete list of simulated parameter values and noise levels is reported in the Supplementary Material.
Interpretation of Positron Emission Tomography Data The PET data are commonly modeled as a convolution of the metabolitecorrected input function and the tracer IRF in the tissue:2 Z t C T ðtÞ ¼ C P IRF ðtÞ ¼ IRFðt - sÞC P ðsÞds ð1Þ 0
where CT(t) is the PET tracer TAC in an ROI or voxel and CP(t) is the metabolite-corrected input function. In the case of 2TCM, IRF is a function of the K1 to k4 values of the corresponding CM, and is indicated here as IRFCM: IRF CM ðtÞ ¼ K 1 Uf ðt; k 2 ; k 3 ; k 4 Þ ¼ K 1 ½πe - tα1 þ ð1 - πÞe - tα2 where α1ð2Þ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 k 2 þ k 3 þ k 4 7 ðk 2 þ k 3 þ k 4 Þ2 - 4k 2 k 4 2
ð2Þ
ð3Þ
and π¼
k 3 þ k 4 - α2 α1 - α2
ð4Þ
However, according to the theory of intravascular indicator dilutions, which was developed in the seminal work of Meier and Zierler22 and then extended to the general case of extravascular tracers (as in the case of PET) by O'Sullivan et al,3 an equation analog to (1) can be written: Z t C T ðtÞ ¼ C P IRFðtÞ ¼ K Rðt - sÞC P ðsÞds ð5Þ 0
Equation (5) provides a probabilistic framework for representation of PET tracer data in terms of a convolution between the metabolite-corrected input function and a tissue residue function R(t). K is a proportionality constant, interpretable as an overall flow. The residue is a scaled survival Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1370 function associated with tracer residence in the tissue. Physically, the residue function describes the fraction of indicator that remains in the tissue in response to an idealized bolus (Dirac delta function) input concentration at time t = 0. Initially, the residue must be unity (R(0) = 1) and from there it decreases (at least does not increase) as t increases. Assuming the residue is differentiable, it can be represented in the form: Z t hðτÞdτ ð6Þ RðtÞ ¼ 1 0
where h(t) is a probability density function. Meier and Zierler22 developed an interpretation of h(t) in terms of the residence time distribution associated with an ensemble of capillary transport paths within the region. Indicator molecules are assumed to randomly sample such paths. Once a path is selected, the travel time of the molecule in the tissue is determined. Either due to flow, volume variations or binding, the transit time of the indicator in each transport path is considered to be distinct. From the distribution of transport paths there arises a distribution of transit times that would be experienced by a randomly chosen indicator molecule. The probability that an indicator molecule ends up on a path with transit time between t and t+Δt is proportional to h(t)Δt as Δt gets small. Hence if a collection of indicator molecules is introduced in trace amounts, so they do not interact, (6) gives the fraction of molecules that experiences transit times greater than t.
Outcome Measures Expressing (1) as in (5), and considering nonparametric inference for the residue, gives rise to a deconvolution problem, providing a novel approach to the analysis of dynamic PET data that is not reliant on specific compartmental modeling assumptions.3 However, once the IRF has been nonparametrically estimated via deconvolution, it remains to establish which are the outcomes that can be actually computed from the IRF, and how they compare with respect to traditional outcomes used in PET. The first outcome measure usually computed in PET for a reversible tracer is VT, whose formulation in terms of the 2TCM K1 to k4 values is:2 Z 1 K1 k3 VT ¼ 1þ f ðt; k 2 ; k 3 ; k 4 Þdt ð7Þ ¼ K1 k2 k4 0 VT equals the integral up to infinity of the corresponding IRF, which is only made possible by the parametric formulation of the IRF. When the IRF is computed nonparametrically, however, the IRF is estimated only up to the end time of the PET scan. A different definition is therefore necessary, as proposed by Jiang et al,4 in which a surrogate ‘volume of distribution’ for the tracer is calculated by integrating the IRF only up to the end time Tend of the scan: Z T end V T - end ¼ K RðπÞdπ ≠V T ð8Þ 0
To obtain a volume that is comparable to the approach used with CMs, O’Sullivan et al3 proposed an alternative formulation for VT based on the idea that, if the residence density is partitioned into quantiles, the volumes occupied by fast, intermediate, and slow moving tracer label can be separately assessed.3 Their formulation uses the scaled first quantile of the residence density: V T - NP1 ¼
KR - 1 ð1=4Þ log ð4Þ
−1
ð9Þ
where R (1/4) indicates the first quantile of the residence, which is the time at which the IRF reaches a fourth of its initial value, and the scaling log(4) ensures that if the residence density is exponential (i.e., as in the 1TCM case), then R − 1(1/4)/log(4) is equal to the mean of the residence distribution. The scaling factor in Equation (9) differs from the one reported in the paper by O’Sullivan et al.3 We closely checked our derivation of the scaling factor under the assumption that the residence density is exponential (h(t) = α × exp(− αt)). We confirmed our derivation using the reported clinical and simulated data. A closed-form expression for the scaling factor when the residence density is biexponential or more cannot be derived. While the use of the first quantile of the residence was suggested for the irreversible tracer [18F]FDG,3 we propose here to further generalize and optimize this nonparametric formulation of VT for reversible tracers, by considering not only the first quantile of the residence density (i.e., R −1(1/4)), but also a large number of quantiles (i.e., R − 1(n), for n = N1, …, NQ), Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
and then averaging over all of them: V T - NP2 ¼
K Xn¼NQ R - 1 ðnÞ n¼N1 log ð1=nÞ N
ð10Þ
This is a natural extension of the idea of tying the nonparametric estimate of the IRF with a particular counterpart; rather than relying only on a single estimated quantile with its inherent uncertainty, this approach will tend to ‘sharpen’ the estimator by considering many estimated quantiles with the expectation that the uncertainties will cancel each other out to some extent. Usually with CMs, once VT is computed according to (7), the binding potential BPP in any ROI/voxel can be indirectly estimated as BPP = VT − VTref, where VTref is the tracer VT in an RR devoid of the target of interest.1 The same formulation of BPP can be adopted for the VT calculated according to Equations (8, 9, and 10), to obtain BPP-end, BPP-NP1, and BPP-NP2, respectively.
Singular Value Decomposition Potentially, any of many nonparametric deconvolution approaches can be used to estimate the IRF starting from known metabolite-corrected input function and tracer TACs in the brain tissue. O’Sullivan et al3 propose a technique based on regularized cubic B-spline approximation of the residence time distribution, but do not claim any optimality of the B-spline estimators. Jiang et al4 use functional principal component analysis. Here, we propose the use of SVD. Specifically, if the tissue concentration CT(t) of the tracer in response to a metabolite-corrected input function CP(t) is given by Equation (5), then assuming that the metabolite-corrected input function and the tracer TACs in the brain tissue are measured at equally spaced time points t1, t2, …, tN, with Δt as the space between time points of measurement, Equation (5) can be discretized, assuming that over short time intervals Δt, the residue function and metabolite-corrected input function are constant in time: Z tj Xj C P ðπÞRðt - πÞdπ ffi KΔt C ðt ÞRðtj - t i Þ ð11Þ C T ðtj Þ ¼ K i¼0 P j 0
or
0
C P ðt1 Þ B C P ðt2 Þ B KΔt@ … C P ðtN Þ
0 C P ðt1 Þ … C P ðtN - 1 Þ
10 1 0 1 Rðt1 Þ C T ðt1 Þ … 0 B Rðt2 Þ C B C T ðt2 Þ C … 0 C CB C¼B C … … A@ … A @ … A RðtN Þ C T ðtN Þ … C P ðt1 Þ
ð12Þ
The requirement for equidistantly sampled data derives from the Shannon’s sampling theorem,23 which implies that any band-limited signal can be numerically reconstructed from its regularly spaced samples. The solution of the nonuniform sampling problem involves much more complex theoretical and computational issues (see Discussion).24 In matrix form, (12) is expressed as: Ab ¼ c
ð13Þ
where b contains the elements of R(ti), i = 1, 2,…, N, and c are the measured tracer TACs in the brain tissue CT(ti), i = 1, 2,…, N. Equation (13) can be solved iteratively for the elements of b. This approach is, however, extremely sensitive to noise, potentially causing the estimated R(t) to oscillate. A technique to solve (13) is SVD, which constructs matrices V, W, and UT, so that the inverse of A in (13) can be written: A - 1 ¼ VWUT
ð14Þ T
where W is a diagonal matrix and V and U are orthogonal matrices. Given this inverse matrix, b, and consequently R(t), is calculated as b ¼ A - 1 c ¼ VWUT c
ð15Þ
The main force of the SVD is that the diagonal elements in W are zero or close to zero corresponding to linear equations in (12) that are close to being linear combinations of each other. This fact thus allows one to identify elements in the matrix A that causes the solution b to oscillate or otherwise be meaningless in a biomedical modeling context. By eliminating (setting equal to zero) diagonal elements below a certain threshold in W, one can consequently minimize these effects before calculating b.25 Here, we propose an automatic data-driven algorithm for selecting the threshold for SVD for each ROI that is based on cross-validation. Specifically, for each TAC in each considered ROI, the threshold is automatically selected by: © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1371 1. For i = 1 to i = N (with N the total number of original sampling time points) (1) Removing the TAC point corresponding to the sampling time point ti; (2) For threshold T = 2% to T = 40% (with step = 1%), applying SVD with threshold value T to the TAC containing the remaining N − 1 data points; (3) Storing the threshold T that gives the best prediction of the TAC data at the time point ti that was left out (i.e., calculated as the minimum of the square difference between original TAC value at the excluded time point ti and the reconvolved TAC obtained with SVD at the same ti); 2. Selecting the average value (across excluded time points ti) of the stored thresholds. The suggested threshold range is a generous range that has proven optimal before for physiologic signals.9,10
Metrics of Comparison for Simulated Data For all simulated cases and noise conditions, we ran SVD with the automatic threshold selection, and calculated VT-NP1, VT-NP2, and VT-end from the deconvolved IRF estimates. BPP-NP1, BPP-NP2, and BPP-end were derived by subtracting VT-NP1, VT-NP2, and VT-end of the RR. We compared the distribution of outcome measures estimates obtained with the true simulated values, and calculated variance (across the 1,000 estimates of the outcome measure in each ROI, the median variance across ROIs is reported) and bias (as distance between median of the 1,000 estimates of the outcome measure in each ROI and the corresponding true value, the median bias across ROIs is then considered) of such outcome measures.
Application to Test–Retest Data For all tracers, subjects, and ROIs, we calculated VT and BPP = VT − VTref using CMs (2TCM for CUMI,14 PE2I,16 and WAY;17 1TCM for DASB15); the RR was the gray matter (for CUMI14 and DASB15), the white matter (for WAY17), or the totality (for PE2I16) of the cerebellum. We also ran SVD with the automatic threshold selection, and calculated VT-NP1, VT-NP2, and VT-end from the deconvolved IRF. BPP-NP1, BPP-NP2, and BPP-end were derived by subtracting VT-NP1, VT-NP2, and VT-end, respectively, of the RR. We first compared VT-NP1, VT-NP2, and VT-end versus VT values using correlation (r), slope and intercept of the regression line, with VT values being the independent variables. BPP-NP1, BPP-NP2, and BPP-end were similarly compared with BPP values. We then compared the test–retest performance of VT, VT-NP1, VT-NP2, VT-end, and BPP, BPP-NP1, BPP-NP2, BPP-end, and using five metrics:15 (1) Percent difference (PD), calculated as the absolute difference between the test and retest values divided by their average (the lowest the PD, the highest the reproducibility); (2) Within-subject of squares (WSMSS): sum P P mean WSMSS ¼ n1 ni¼1 2j¼1 X ij - X i , where Xi1 represents the ‘test’ value of an outcome measure for subject i, Xi2 the ‘retest’ value for the same subject, n is the number of pairs of subjects, and X i ¼ ðX i1 þ X i2 Þ=2. This represents the variance in the estimation of a subject-specific outcome measure (the lowest the WSMSS, the highest the reproducibility); BSMSS - WSMSS (3) Interclass correlation coefficient (ICC): ICC ¼ BSMSSþðk - 1ÞWSMSS, where BSMSS is the between-subject mean sum of squares
P 1 Pk and k is the number of repeated X X BSMSS ¼ n1 2n ij i i¼1 k j¼1 observations (k = 2 in this study; − 1 corresponds to no reliability, 1 to maximum reliability); (4) Variance of an outcome measure across subjects (sample variance, across subjects, of all observations for each ROI); and (5) Identifiability based on bootstrap resampling of residuals; we computed 100 bootstrap samples26 of the TAC data for each subject’s ROI, estimated the outcome measures for each of these samples, and measured the variability of these estimates using the robust median absolute deviation criterion: median [|Zj – median(Z1, Z2, y , Z100)|, j = 1,2, y ,100]. © 2015 ISCBFM
RESULTS Test–Retest Data Comparison with CMs estimates. Figures 1 and 2 show for all tracers, subjects, and ROIs, the scatter plots of VT-NP1, VT-NP2, and VT-end, and BPP-NP1, BPP-NP2, and BPP-end, respectively, obtained using SVD versus the VT and BPP values obtained using CMs. The identity and regression lines, and the numerical results of the regression analysis, are reported in each plot. Correlation is high (40.90) in all cases, with exception of BPP-NP1 for PE2I (r = 0.740) and BPP-end for CUMI (r = 0.875). However, the slope values differ across outcome measures and tracers. As expected from its formulation in comparison with VT, VT-end underestimation increases as the value of the estimated outcome increases. VT-end shows the lowest slope values across the different outcomes for all tracers: the minimum value is 0.466 (PE2I) and the maximum is 0.619 (DASB). VT-NP1 and VT-NP2 show similar results with each other, and the slope values are closer to unity than with VT-end (VT-NP1: minimum value is 0.511, PE2I, maximum is 1.019, WAY; VT-NP2: minimum value is 0.827, PE2I, maximum is 0.944, DASB). As expected by definition, VT-NP2 shows higher correlation than VT-NP1 for all tracers, while the slope is closer to unity than with VT-NP1 for DASB and PE2I. For all tracers, slope values for BPP decrease with respect to the comparison with VT values. In general, BPP-NP2 shows the best correlation values for all tracers (0.922 for CUMI, 0.986 for DASB, 0.973 for PE2I, and 0.942 for WAY), and the best slope values for DASB (0.928) and PE2I (0.794). BPP-NP1 shows the best slope value for CUMI (0.858) and WAY (1.044). Test–retest metrics. Figures 3 and 4 show the boxplots of the five considered metrics for all tracers and formulations of the outcome measure (VT, VT-NP1, VT-NP2, VT-end and BPP, BPP-NP1, BPP-NP2, BPP-end, respectively). In each box, the central mark is the median value of the metric (across subjects/pairs and ROIs), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. A complete list of the numerical results corresponding to Figures 3 and 4 is reported in the Supplementary Material. Table 1 summarizes which outcome measure provides the best performance. For the clinical data, the outcome measures with the best correlation and slope (with the corresponding outcome measures estimated with CMs) are reported separately from the test–retest metrics; the best performance in the test–retest metrics is reported both considering and not considering VT-end (or BPP-end). VT-end shows the lowest WSMSS and variance, and the best identifiability, for all tracers. VT-end also shows the lowest PD for CUMI, DASB, and PE2I, and the highest ICC for CUMI and DASB. Correspondingly, BPP-end shows the lowest WSMSS and variance for all tracers. BPP-end also shows the highest ICC for all tracers, besides WAY, the lowest PD for DASB and PE2I, and the best identifiability for all tracers, with exclusion of DASB. VT and BPP show the lowest PD and the highest ICC for WAY. BPP also shows the best identifiability for DASB. VT-NP1 and BPP-NP1 show the worst identifiability for all tracers. BPP-NP1 also shows the lowest ICC for all tracers, besides WAY. VT-NP2 shows the highest ICC for PE2I (0.80), while BPP-NP2 shows the lowest PD for CUMI. Simulated Data Figures 5 and 6 show the distribution of outcome measures obtained using SVD in the 1TCM and 2TCM case, respectively (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2), with simulated noise level corresponding to α = 1. Each region is Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1372
Figure 1. Scatter plots of VT-NP1, VT-NP2, and VT-end obtained by singular value decomposition (SVD) versus the VT values obtained using compartment models (CMs) for all considered tracers, subjects, and regions of interest (ROIs). The identity (solid) and regression (dotted) lines are reported in each plot, together with the numerical results of the regression analysis. VT, total volume of distribution.
reported separately (the cerebellum gray matter is the RR for calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation of the outcome measure, while the black solid vertical line represents the true simulated value for that outcome measure and region. In the 1TCM case, VT-end always underestimates VT the most across the three formulations, except for the cerebellum gray matter, although with much less variability. In the 2TCM case, VT-end performs the best for the three regions with the lowest VT (cerebellum gray matter, occipital lobe, and cingulate), but underestimates the most the region with the highest VT (hippocampus). Correspondingly, BPP-end always underestimates BPP the most in both cases. In the 1TCM case, VT-NP1 and VT-NP2 show similar median values for VT, with VT-NP2 performing better than VT-NP1 for the three Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
regions with highest simulated VT (dorsal caudate, amydgala, ventral striatum). In the 2TCM case, it is VT-NP1 that slightly outperforms VT-NP2 in all regions in terms of bias, but not variance of the estimates distribution. Corresponding BPP-NP1 and BPP-NP2 values are similar in the 1TCM case, while BPP-NP1 tends to overestimate the BPP in all the regions in the 2TCM case. BPP-NP2 shows the smallest bias across regions and less dispersed distributions than the ones of BPP-NP1 in the 2TCM case. As reported in Table 1, with the noise present at the ROI level (α = 1), VT-end and BPP-end show the best variance in both simulated cases: 0.039 and 0.084 (1TCM); 0.043 and 0.090 (2TCM). VT-end shows the smallest bias (−0.298) for 2TCM, while VT-NP2 shows the smallest bias for 1TCM (−0.426), and BPP-NP2 shows the smallest bias for both 1TCM (0.129) and 2TCM (0.351). The results relative to the noise analysis are reported in the Supplementary Material. In both 1TCM and 2TCM cases, the bias of © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1373
Figure 2. Scatter plots of BPP-NP1, BPP-NP2, and BPP-end obtained by singular value decomposition (SVD) versus the BPP values obtained using compartment models (CMs) for all considered tracers, subjects, and regions of interest (ROIs). The identity (solid) and regression (dotted) lines are reported in each plot, together with the numerical results of the regression analysis. BPp, binding potential.
the VT-end and BPP-end estimates is almost constant across the five noise conditions. The corresponding variance, although it increases with the noise in the data, remains significantly lower than the corresponding variance of the other two nonparametric formulations for the same outcome measure. In both 1TCM and 2TCM cases, the pattern with which the bias of VT-NP1 and BPP-NP1, and VT-NP2 and BPP-NP2, changes with noise is similar, with both formulations that tend to overestimate the true value at elevated noise levels. In the 1TCM case, the bias of VT-NP2 and BPP-NP2 increases the most with the noise level. The corresponding variance values increase with noise, with VT-NP1 and BPP-NP1 providing the worst performance in the 2TCM case. In the 1TCM case, the variance of VT-NP1 and VT-NP2, and BPP-NP1 and BPP-NP2, is almost identical. © 2015 ISCBFM
DISCUSSION We applied nonparametric deconvolution to estimate the IRF of PET reversible tracers in several regions, using both simulated and clinical test–retest data sets, computed several functionals of the IRF to estimate tracer VT and BPs, and extensively compared these functionals with traditional outcome measures derived using CMs. We showed that nonparametric deconvolution, which does not require any model assumptions, and the functionals that can be computed from the recovered IRF can provide estimates of VT and BPP that are not only very close to, and highly correlated with, the estimated obtained with traditional CMs, but can also show better test–retest performance than outcome measures derived from CMs. Therefore, nonparametric deconvolution and the proposed functionals offer a promising alternative to CMs for the Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1374
Figure 3. Boxplots of the five considered test–retest metrics for all considered tracers and the four different estimations of the outcome measure (VT, VT-NP1, VT-NP2, VT-end). In each box, the central mark is the median value of the metric (across subjects/pairs and regions), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. ICC, interclass correlation coefficient; PD, percent difference; VT, total volume of distribution; WSMSS, within-subject mean sum of squares.
quantification of PET dynamic data, especially when CMs assumptions do not hold or when CMs provide poor fits and/or unphysiologic estimates of the outcome measures for a specific PET tracer. The optimized formulation of VT we proposed in Equation (10), and the corresponding BPP-NP2 that can be derived, provide the best performance in terms of bias with respect to, and correlation with, outcome measures with traditional CMs in the majority of the simulated and clinical considered data sets (Table 1). VT-end as in Equation (8), and the corresponding BPP-end that can be derived, provide the best test–retest performance across metrics and, overall, considered tracers; outcome measures derived with traditional CMs show superior test–retest performance than VT-end and BPP-end only for some of the metrics in the case of DASB Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
(identifiability of BPP), PE2I (ICC of VT), and WAY (PD and ICC of VT and BPP). Therefore, VT-end and BPP-end are stable (with respect to noise in the TACs) and reliable outcome measures that could be used, for example, to discriminate different clinical populations from each other and from healthy controls. However, as one would expect from their definitions, results on simulated and clinical data sets suggest that VT-end and BPP-end often provide the greatest underestimation of VT and BPP (with the exception of PE2I BPP, for which BPP-NP1 shows the greatest underestimation). This is perhaps less of an issue if the purpose of the study is to examine differences between groups (or across values of some continuous variable) but if the primary goal of the analysis is estimation of standard outcome measures, the bias of these estimates must be considered. In such a situation, VT-NP1 and/or VT-NP2 (and the © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1375
Figure 4. Boxplots of the five considered test–retest metrics for all considered tracers and the four different estimations of the outcome measure (BPP, BPP-NP1, BPP-NP2, BPP-end). In each box, the central mark is the median value of the metric (across subjects/pairs and regions), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. BP, binding potential; ICC, interclass correlation coefficient; PD, percent difference; WSMSS, within-subject mean sum of squares.
corresponding BPP-NP1 and BPP-NP2) could be considered, as they show better test–retest performance than CMs outcomes according to the majority of the considered metrics for all tracers, with the exception of WAY. An intrinsic limitation of the nonparametric approach is that the IRF can be estimated only up to the end of the PET scan, and no inference up to infinity can therefore be made without some strong assumptions. For example, it may be the case that, for tracers that slowly reach equilibrium, the nonparametric estimation of VT and BPP according to Equations (9 and 10) is biased by the fact that, by the end of the scan, the IRF has not yet reached a certain fraction of its initial value, and therefore the corresponding kth quantile of the residence cannot be determined. Estimation of VT-end and BPP-end instead might be preferred in these cases, with proper interpretation of the results. © 2015 ISCBFM
We have implemented here nonparametric deconvolution using SVD9,10 with an automatic data-driven selection of the threshold for each ROI. This is a fast, easy to implement, and widely adopted approach to deconvolution that easily lends itself to implementation at the voxel level. However, we claim no optimality in either the choice of the deconvolution approach or the selection of the thresholds. It is known that SVD performance can be sensitive to the level of the noise present in the data.27,28 Our simulation results at different noise level (see Supplementary material) suggest that, at the voxel level, or for particularly noisy tracers, either VT-end and BPP-end should be preferred as outcome measures because they are less affected by noise, or more robust approaches to deconvolution, such as nonlinear stochastic regularization28 or Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1376 Table 1.
Summary of which outcome measure provides the best performance in each of the simulated and clinical data sets Simulated data
1TC 2TC Clinical data CUMI DASB PE2I WAY
VT Best variance
Best bias
BPp Best variance
Best bias
VT-end VT-end
VT-NP2 VT-end
BPp-end BPp-end
BPp-NP2 BPp-NP2
Best correlation
Best slope
Best correlation
Best slope
VT-NP2 VT-NP2 VT-NP2 VT-NP2
VT-NP1 VT-NP2 VT-NP2 VT-NP2
BPp-NP2 BPp-NP2 BPp-NP2 BPp-NP2
BPp-NP1 BPp-NP2 BPp-NP2 BPp-NP2
Test–retest metrics
VT (considering VT-end) CUMI DASB PE2I WAY VT (not considering VT-end) CUMI DASB PE2I WAY BPp (considering BPp-end) CUMI DASB PE2I WAY BPp (not considering BPp-end) CUMI DASB PE2I WAY
PD
WSMSS
ICC
Variance
Identifiability
VT-end VT-end VT-end VT
VT-end VT-end VT-end VT-end
VT-end VT-end VT/VT-NP2 VT
VT-end VT-end VT-end VT-end
VT-end VT-end VT-end VT-end
VT-NP2 VT VT-NP1 VT
VT-NP2 VT-NP2 VT-NP1 VT
VT-NP2 VT/VT-NP2 VT/VT-NP2 VT
VT-NP2 VT-NP1 VT-NP1 VT
VT-NP2 VT VT-NP2 VT
BPp-NP2 BPp-end BPp-end BPp
BPp-end BPp-end BPp-end BPp-end
BPp-end BPp-end BPp-end BPp
BPp-end BPp-end BPp-end BPp-end
BPp-end BPp BPp-end BPp-end
BPp-NP2 BPp BPp BPp
BPp-NP2 BPp-NP2 BPp-NP1 BPp
BPp-NP2 BPp-NP2 BPp-NP2 BPp
BPp-NP2 BPp BPp-NP1 BPp
BPp BPp BPp-NP2 BPp
Abbreviations: BPp, binding potential; CM, compartment model; ICC, interclass correlation coefficient; PD, percent difference; 1TC, one-tissue compartment; 2TC, two-tissue compartment; VT, total volume of distribution; WSMSS, within-subject mean sum of squares. For the simulated data, the outcome measures with the best variance and best bias (with respect to the true value of the outcome measure) are reported. For the clinical data, the outcome measures with the best correlation and slope (with the corresponding outcome measures estimated with CM) are reported separately from the test–retest metrics; the best performance in the test–retest metrics is reported considering and not considering VT-end (or BPP-end).
kernel-based linear system identification,29,30 should be considered for reconstructing the IRF. The cross-validation algorithm we propose here for thresholding the singular values in the decomposition to constrain the propagation of noise from the data to the IRF estimates is data driven and has proven to be flexible enough to provide reasonably regularized solutions for the considered four tracers and simulated data. However, other strategies exist,7,31–33 which may provide improved solutions. The optimization of the thresholding algorithm, which likely depends on the tracer at hand, is however beyond the scope of this article. Singular value decomposition, as many other nonparametric deconvolution approaches, requires the input data to be equidistantly sampled.23 The PET data involve radioactive decay, and are usually acquired using increasing frame lengths to ensure a comparable number of counts and thus statistical quality per data point. This implies that the PET data need to be interpolated on an equally-spaced grid of temporal values before applying SVD. The possible effect of the interpolation on the estimation of VT and BPP should be further investigated, especially when using SVD at the voxel level where, due to the large number of voxels to quantify, adopting a very fine interpolation grid can Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
lead to out-of-memory issues during the quantification process. While Shannon’s sampling theorem23 implies that any band-limited signal can be reconstructed from its regularly spaced samples, the solution of the nonuniform sampling problem poses more difficulties.24 There exist a vast literature on nonuniform sampling theory34,35 and methods claiming to efficiently reconstruct a function from its nonuniform samples,36–38 including nonuniform SVD.39 The investigation of the performance on PET data of nonuniform deconvolution approaches is however beyond the scope of this article. It is important to state that comparing the performance of the proposed nonparametric outcome measures with the corresponding CMs estimates is not generally useful when the CM is not valid for the tracer in question. Here, we have considered four tracers for which we have already assessed the validity of the CMs,14–17 so the proposed comparison is meaningful, although it may be questioned in certain individual regions. Certainly, the possibility of quantifying the binding of PET tracers model-free and nonparametrically is particularly attractive for those tracers for which the assumptions of CMs do not hold or CMs do not provide a good description of the data. In these cases, however, the CMs estimates of the binding cannot be used as a © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1377
Figure 5. Distribution of outcome measures estimated using singular value decomposition (SVD) and the three nonparametric formulations (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2) across the 1,000 simulated one-tissue compartment model (1TCM) cases in each region; the six simulated regions are reported separately (the cerebellum gray matter is used as the reference region (RR) for the calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation, while the black solid vertical line represents the true simulated value for that outcome measure and region. The data underlying the simulations were derived from clinical [11C]DASB scans modeled with 1TCM. BPp, binding potential; VT, total volume of distribution.
reference to evaluate the performance of the nonparametric approach, and validation with in vitro results would be recommended.
© 2015 ISCBFM
Future directions of investigation also include the extension of the presented nonparametric approach to quantification in the absence of a metabolite-corrected input function.
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1378
Figure 6. Distribution of outcome measures estimated using singular value decomposition (SVD) and the three nonparametric formulations (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2) across the 1,000 simulated two-tissue compartment model (2TCM) cases in each region; the five simulated regions are reported separately (the cerebellum gray matter is used as the reference region (RR) for the calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation, while the black solid vertical line represents the true simulated value for that outcome measure and region. The data underlying the simulations were derived from clinical [11C]CUMI scans modeled with 2TCM. BPp, binding potential; VT, total volume of distribution.
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
© 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1379 CONCLUSIONS Nonparametric deconvolution does not require any model assumptions and allows for estimates of VT and BPP that are very close to the estimates obtained with traditional CMs, and can also show superior test–retest performance. Nonparametric deconvolution and the proposed functionals therefore offer a viable alternative to CMs for the quantification of PET dynamic data, especially when CMs assumptions do not hold or CMs provide poor fits and/or unphysiologic estimates of the outcome measures for a specific PET tracer. AUTHOR CONTRIBUTIONS Dr Zanderigo was responsible for ideation of the methodologies, drafting of the manuscript, creation of the software, analysis and interpretation of the data. Dr Parsey contributed to conception and revision of the manuscript, and interpretation of the data. Dr Ogden contributed to ideation of the methodologies, interpretation of the data, critical revision of the manuscript, and final approval.
DISCLOSURE/CONFLICT OF INTEREST The authors declare no conflict of interest.
REFERENCES 1 Innis RB, Cunningham VJ, Delforge J, Fujita M, Gjedde A, Gunn RN et al. Consensus nomenclature for in vivo imaging of reversibly binding radioligands. J Cereb Blood Flow Metab2007; 27: 1533–1539. 2 Slifstein M, Laruelle M. Models and methods for derivation of in vivo neuroreceptor parameters with PET and SPECT reversible radiotracers. Nucl Med Biol2001; 28: 595–608. 3 O'Sullivan F, Muzi M, Spence AM, Mankoff DM, O'Sullivan JN, Fitzgerald N et al. Nonparametric residue analysis of dynamic PET data with application to cerebral FDG studies in normals. J Am Stat Assoc 2009; 104: 556–571. 4 Jiang C-R, Aston JAD, Wang J-L. A functional approach to deconvolve dynamic neuroimaging data 2014; arXiv.org/stat/arXiv:1411.2051. 5 Cunningham VJ, Jones T. Spectral analysis of dynamic PET studies. J Cereb Blood Flow Metab 1993; 13: 15–23. 6 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: 1425–1439. 7 Turkheimer FE, Hinz R, Gunn RN, Aston JA, Gunn SR, Cunningham VJ. Rankshaping regularization of exponential spectral analysis for application to functional parametric mapping. Phys Med Biol 2003; 48: 3819–3841. 8 Schmidt K. Which linear compartmental systems can be analyzed by spectral analysis of PET output data summed over all compartments? J Cereb Blood Flow Metab 1999; 19: 560–569. 9 Ostergaard L, Weisskoff RM, Chesler DA, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis. Magn Reson Med 1996; 36: 715–725. 10 Ostergaard L, Sorensen AG, Kwong KK, Weisskoff RM, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part II: Experimental comparison and preliminary results. Magn Reson Med 1996; 36: 726–736. 11 Buchert R, Wilke F, van den Hoff J, Mester J. Improved statistical power of the multilinear reference tissue approach to the quantification of neuroreceptor ligand binding by regularization. J Cereb Blood Flow Metab 2003; 23: 612–620. 12 Constantinescu CC, Bouman C, Morris ED. Nonparametric extraction of transient changes in neurotransmitter concentration from dynamic PET data. IEEE Trans Med Imaging 2007; 26: 359–373. 13 Grecchi E, Doyle OM, Bertoldo A, Pavese N, Turkheimer FE. Brain shaving: adaptive detection for brain PET data. Phys Med Biol 2014; 59: 2517–2534. 14 Milak MS, DeLorenzo C, Zanderigo F, Prabhakaran J, Kumar JS, Majo VJ et al. In vivo quantification of human serotonin 1A receptor using 11C-CUMI-101, an agonist PET radiotracer. J Nucl Med 2010; 51: 1892–1900.
15 Ogden RT, Ojha A, Erlandsson K, Oquendo MA, Mann JJ, Parsey RV. In vivo quantification of serotonin transporters using [(11)C]DASB and positron emission tomography in humans: modeling considerations. J Cereb Blood Flow Metab 2007; 27: 205–217. 16 DeLorenzo C, Kumar JS, Zanderigo F, Mann JJ, Parsey RV. Modeling considerations for in vivo quantification of the dopamine transporter using [(11)C]PE2I and positron emission tomography. J Cereb Blood Flow Metab 2009; 29: 1332–1345. 17 Parsey RV, Slifstein M, Hwang DR, Abi-Dargham A, Simpson N, Mawlawi O et al. Validation and reproducibility of measurement of 5-HT1A receptor parameters with [carbonyl-11C]WAY-100635 in humans: comparison of arterial and reference tisssue input functions. J Cereb Blood Flow Metab 2000; 20: 1111–1133. 18 Jenkinson M, Smith S. A global optimisation method for robust affine registration of brain images. Med Image Anal 2001; 5: 143–156. 19 Smith SM. Fast robust automated brain extraction. Hum Brain Mapp 2002; 17: 143–155. 20 Parsey RV, Ogden RT, Mann JJ. Determination of volume of distribution using likelihood estimation in graphical analysis: elimination of estimation bias. J Cereb Blood Flow Metab 2003; 23: 1471–1478. 21 Parsey RV, Ojha A, Ogden RT, Erlandsson K, Kumar D, Landgrebe M et al. Metabolite considerations in the in vivo quantification of serotonin transporters using 11C-DASB and PET in humans. J Nucl Med 2006; 47: 1796–1802. 22 Meier P, Zierler KL. On the theory of the indicator-dilution method for measurement of blood flow and volume. J Appl Physiol 1954; 6: 731–744. 23 Shannon C. A mathematical theory of communication. Bell Syst Tech J 1948; 27: 379–623. 24 Strohmer T. Numerical analysis of the non-uniform sampling problem. J Comput Appl Math 2000; 122: 297–316. 25 Press WH, Teukolsky SA, Vetterling WT, Flannery BT. Numerical Recipes in C. 2nd edition Cambridge University Press: Oxford, 1992. 26 Ogden RT, Tarpey T. Estimation in regression models with externally estimated parameters. Biostatistics 2006; 7: 115–129. 27 Knutsson L, Stahlberg F, Wirestam R. Aspects on the accuracy of cerebral perfusion parameters obtained by dynamic susceptibility contrast MRI: a simulation study. Magn Reson Imaging 2004; 22: 789–798. 28 Zanderigo F, Bertoldo A, Pillonetto G, Cobelli Ast C. Nonlinear stochastic regularization to characterize tissue residue function in bolus-tracking MRI: assessment and comparison with SVD, block-circulant SVD, and Tikhonov. IEEE Trans Biomed Eng 2009; 56: 1287–1297. 29 Pillonetto G. DNG. A new kernel-based approach for linear system identification. Automatica 2010; 46: 81–93. 30 Castellaro M, Peruzzo D, Mehndiratta A, Pillonetto G, Petersen ET, Golay X et al. Estimation of arterial arrival time and cerebral blood flow from QUASAR arterial spin labeling using stable spline. Magn Reson Med 2014; doi: 10.1002/mrm.25525 (e-pub ahead of print). 31 Murase K, Yamazaki Y, Miyazaki S. Deconvolution analysis of dynamic contrastenhanced data based on singular value decomposition optimized by generalized cross validation. Magn Reson Med Sci2004; 3: 165–175. 32 Sourbron S, Luypaert R, Van Schuerbeek P, Dujardin M, Stadnik T. Choice of the regularization parameter for perfusion quantification with MRI. Phys Med Biol2004; 49: 3307–3324. 33 Liu HL, Pu Y, Liu Y, Nickerson L, Andrews T, Fox PT et al. Cerebral blood flow measurement by dynamic contrast MRI using singular value decomposition with an adaptive threshold. Magn Reson Med 1999; 42: 167–172. 34 Landau HJ. Necessary density conditions for sampling and interpolation of certain entire functions. Acta Math 1967; 117: 37–52. 35 Feichtinger HG, Gröchenig KH. Theory and practice of irregular sampling. Benedetto J, Frazier M, (Eds.), Wavelets: Mathematics and Applications. CRC Press: Boca Raton, FL, 1994; 305–363. 36 Yao K TJ. On some stability and interpolatory properties of nonuniform sampling expansions. IEEE Trans Circuit Theory 1967; 14: 404–408. 37 Feichtinger HG GK, Strohmer T. Efficient numerical methods in non-uniform sampling theory. Numer Math 1995; 69: 423–440. 38 De Nicolao G, Liberati D, Sartorio A. Deconvolution of infrequently sampled data for the estimation of growth hormone secretion. IEEE Trans Biomed Eng 1995; 42: 678–687. 39 Wingham D. The reconstruction of a band-limited function and its Fourier transform from a finite number of samples at arbitrary locations by singular value decomposition. IEEE Trans Circuit Theory 1967; 40: 559–570.
Supplementary Information accompanies the paper on the Journal of Cerebral Blood Flow & Metabolism website (http://www.nature. com/jcbfm)
© 2015 ISCBFM
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
ORIGINAL ARTICLE
Model-free quantification of dynamic PET data using nonparametric deconvolution Francesca Zanderigo1, Ramin V Parsey2 and R Todd Ogden3,4,5 Dynamic positron emission tomography (PET) data are usually quantified using compartment models (CMs) or derived graphical approaches. Often, however, CMs either do not properly describe the tracer kinetics, or are not identifiable, leading to nonphysiologic estimates of the tracer binding. The PET data are modeled as the convolution of the metabolite-corrected input function and the tracer impulse response function (IRF) in the tissue. Using nonparametric deconvolution methods, it is possible to obtain model-free estimates of the IRF, from which functionals related to tracer volume of distribution and binding may be computed, but this approach has rarely been applied in PET. Here, we apply nonparametric deconvolution using singular value decomposition to simulated and test–retest clinical PET data with four reversible tracers well characterized by CMs ([11C]CUMI-101, [11C]DASB, [11C]PE2I, and [11C]WAY-100635), and systematically compare reproducibility, reliability, and identifiability of various IRF-derived functionals with that of traditional CMs outcomes. Results show that nonparametric deconvolution, completely free of any model assumptions, allows for estimates of tracer volume of distribution and binding that are very close to the estimates obtained with CMs and, in some cases, show better test–retest performance than CMs outcomes. Journal of Cerebral Blood Flow & Metabolism (2015) 35, 1368–1379; doi:10.1038/jcbfm.2015.65; published online 15 April 2015 Keywords: binding; full quantification; model-free; PET; singular value decomposition
INTRODUCTION Positron emission tomography (PET) is a nuclear imaging technology that uses radioactively labeled molecules (tracers) that are injected into the body and bind, although not exclusively, to a specific target to quantify in vivo biological processes, such as blood flow, metabolism, and the distribution of proteins in the brain and the body. PET is an invaluable tool to establish the genetic and cellular basis of diseases. Fully quantitative PET allows estimation of outcome measures such as, in the case of a reversible tracer, the total volume of distribution (VT) of the injected tracer and a series of binding potentials (BPs) between tracer and target, to indirectly determine the available amount of the target.1 These BPs are in various ways related to the density of available target (Bmax) through the affinity of the tracer for the available binding sites (1/KD).1 Quantification of these outcome measures relies on describing the PET signal associated with the tracer in the tissue (i.e., the time activity curve, TAC) as the convolution between the concentration of free (not protein bound) unmetabolized tracer in the arterial plasma (i.e., the metabolite-corrected input function) and a function specific to the tissue, the impulse response function (IRF).2 The metabolite-corrected input function represents how much of the injected tracer is actually available to bind with the target. IRF is estimated from measures of the input and the tissue PET signal, and would not otherwise be known. Traditionally, PET quantification methods have used parametric approaches, which require specifying a certain parametric form for the IRF (e.g., a mathematical function derived from a physiologic
model of the system, or that approximates the system behavior). This approach therefore relies on fairly strong assumptions, which are not always substantiated in reality, and a priori knowledge of the system under investigation. Traditional PET parametric analysis quantifies these outcome measures by interpreting the PET signal in the brain tissue using compartment models (CMs) to describe the system IRF (the kinetics) of the tracer, in particular a reversible or irreversible two-tissue compartment model (2TCM).1 The tracer signal in the brain tissue is modeled as the sum of two components: part of the tracer is specifically bound to the target of interest (specific binding); part is nonspecifically bound to macromolecules other than the target, or free in the tissue water (nonspecific binding).1,2 Identification of this model involves estimating four rate constants (K1, k2, k3, k4), which describe how quickly the tracer transits from ‘being in the plasma’, to ‘being free to bind in the tissue’, to ‘being specifically bound to the target’.1,2 However, these constants are not always all identifiable, and as a result directly estimating the specific binding component in the model can often be unstable. Therefore, full quantification of this model also requires identifying a reference region (RR) in the brain that is completely devoid of the target.1,2 In a proper (i.e., measurable, valid, reliable, independent of treatment effects and groups) RR devoid of the target of interest, the tracer can only be nonspecifically bound to macromolecules other than the target,1,2 and its kinetics can be described with a one-tissue compartment model1 (1TCM) regulated by only two rate constants (K1’ and k2’). The tracer volume of distribution in such a region is assumed to represent (and used to estimate) a common nonspecific binding
1 Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute and Columbia University, New York, New York, USA; 2Department of Psychiatry and Radiology, Stony Brook University, Stony Brook, New York, USA; 3Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute, New York, New York, USA; 4Department of Psychiatry, Columbia University, College of Physicians and Surgeons, New York, New York, USA and 5Department of Biostatistics, Columbia University, Mailman School of Public Health, New York, New York, USA. Correspondence: Dr F Zanderigo, Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute and Columbia University, 1051 Riverside Drive, New York, NY 10032, USA. E-mail: [email protected] or [email protected] Received 1 December 2014; revised 25 February 2015; accepted 18 March 2015; published online 15 April 2015
Model-free nonparametric PET quantification F Zanderigo et al
1369 across the brain, and helps isolate the component that is only specifically bound to the target. Often, however, CMs or graphical approaches derived from CMs2 either do not properly describe the tracer kinetics, or are not identifiable, leading to nonphysiologic estimates of the VT and of the binding to the target. On the other hand, nonparametric deconvolution approaches do not impose any specific formulation or model for the IRF, and only rely on very loose requirements, often matched in the case of physiologic systems, for the IRF (e.g., degree of smoothness; monotonicity). Nonparametric deconvolution approaches, therefore, allow obtaining model-free estimates of the IRF, from which functionals related to tracer VT and BPs can be computed. Nevertheless, they have been rarely adopted in PET. O'Sullivan et al3 proposed the use of nonparametric residue analysis and demonstrated it with simulations and data using the irreversible tracer [18F]FDG. Only recently, Jiang et al4 proposed a nonparametric approach based on functional principal component analysis and demonstrate it with simulations and data with the tracer [11C]-diprenorphine. Several approaches have been developed that fall between these two extremes.5–7 Generally, these involve expressing the IRF in terms of basis functions that, although derived from the principles of linear time-invariant systems and based on CMs,8 do not impose any predefined compartmental structure. Similar to what we propose using nonparametric deconvolution approaches, these basis function approaches allow for calculation of functionals that can be computed from the IRF as PET outcome measures. Here, we propose the use of nonparametric deconvolution via singular value decomposition (SVD)9,10 to compute functionals of the IRF of PET reversible tracers that are associated with tracer VT and BPs, and extensively compare accuracy, reproducibility, reliability, and identifiability of these different functionals with traditional outcome measures derived using CMs. SVD is an easy-to-implement nonparametric deconvolution approach that is widely used, for example, in the quantification of Dynamic Susceptibility Contrast-Magnetic Resonance Imaging data9,10 and has been used before to obtain regularized least-square solutions.7,11–13 We show relative advantages and disadvantages of using outcome measures derived from either nonparametric deconvolution or traditional parametric approaches, using simulated data and four test–retest clinical datasets with the tracers [11C]CUMI-101 (CUMI), [11C]DASB (DASB), [11C]PE2I (PE2I), and [11C] WAY-100635 (WAY). MATERIALS AND METHODS Subjects Four previously acquired test-retest datasets were considered: 12 scans with CUMI,14 20 with DASB,15 14 with PE2I,16 and 10 with WAY.17 The studies were performed in accordance with the Declaration of Helsinki and The Institutional Review Boards of Columbia University Medical Center and New York State Psychiatric Institute approved the protocol. Subjects gave written informed consent after receiving an explanation of the study.
Positron Emission Tomography Protocol Tracers preparation, emission data acquisition and reconstruction, and determination and fit of the metabolite-corrected input function were obtained as previously described.14–17
Image Analysis Images were analyzed using Matlab Release 2013b (The Mathworks, MA, USA) with extensions to FLIRT version 5.2 (ref. 18) and BET version 1.2 (ref. 19). Motion correction was applied. An RR and several anatomic target regions of interest (ROIs) were traced based on brain atlases and published reports:20,21 13 ROIs for CUMI (hippocampus, entorhinal cortex, insula, posterior parahippocampal gyrus, temporal lobe, amygdala, medial © 2015 ISCBFM
prefrontal cortex, cingulate, orbital, dorsolateral prefrontal cortex, parietal, raphe, occipital lobe; RR: cerebellar gray matter); 6 ROIs for DASB (midbrain, anterior cingulate, amygdala, dorsal putamen, hippocampus, thalamus; RR: cerebellar gray matter); 4 ROIs for PE2I (dorsal caudate, dorsal putamen, ventral striatum; RR: cerebellum); and 13 ROIs for WAY (raphe; anterior cingulate; amygdala; cingulate, dorsolateral prefrontal cortex, hippocampus, insula, medial prefrontal cortex, occipital lobe, orbital, parietal, posterior parahippocampal gyrus, temporal lobe; RR: cerebellar white matter). Selecting these ROIs restricted the analysis to regions for which the binding is expected to be relatively high (with the exception of the RR).
Simulated Data In addition to the clinical test–retest data sets, we considered simulations designed to generate data that imitate characteristics of real data in the case of 1TCM and 2TCM, respectively. For the 1TCM, we used the metabolite-corrected input function measured in one of the considered DASB scans, a tracer whose kinetics is well described by a 1TCM,15 and set the K1 and k2 values in six selected regions (cerebellar gray matter, temporal lobe, hippocampus, dorsal caudate, amygdala, and ventral striatum) based on existing data on DASB healthy control subjects.15 For the 2TCM, we used the metabolitecorrected input function measured in one of the considered CUMI scans, a tracer whose kinetics is well described by a 2TCM,14 and set the K1 to k4 values in five selected regions (cerebellar gray matter, hippocampus, temporal lobe, occipital lobe, and cingulate) based on existing data on CUMI healthy control subjects.14 In both cases, we simulated Gaussian noise with zero mean and standard deviation (SD) calculated as SD(ti) = α. σ/w(ti), where w(ti) is the root square of the duration of the time frame during which the PET data corresponding to time ti has been acquired, and σ indicates the level of noise that we derived from the weighted residual sum of square values obtained by fitting DASB and CUMI scans, respectively, with CMs. Increasing α allows us to study the various estimation techniques in situations with proportionately more noise than what we encounter in our real-data situations. We simulated five noise conditions setting α equal to 1, 2.5, 5, 7.5, and 10, respectively. In all cases and noise conditions, we simulated 1,000 TACs for each of the selected regions. The complete list of simulated parameter values and noise levels is reported in the Supplementary Material.
Interpretation of Positron Emission Tomography Data The PET data are commonly modeled as a convolution of the metabolitecorrected input function and the tracer IRF in the tissue:2 Z t C T ðtÞ ¼ C P IRF ðtÞ ¼ IRFðt - sÞC P ðsÞds ð1Þ 0
where CT(t) is the PET tracer TAC in an ROI or voxel and CP(t) is the metabolite-corrected input function. In the case of 2TCM, IRF is a function of the K1 to k4 values of the corresponding CM, and is indicated here as IRFCM: IRF CM ðtÞ ¼ K 1 Uf ðt; k 2 ; k 3 ; k 4 Þ ¼ K 1 ½πe - tα1 þ ð1 - πÞe - tα2 where α1ð2Þ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 k 2 þ k 3 þ k 4 7 ðk 2 þ k 3 þ k 4 Þ2 - 4k 2 k 4 2
ð2Þ
ð3Þ
and π¼
k 3 þ k 4 - α2 α1 - α2
ð4Þ
However, according to the theory of intravascular indicator dilutions, which was developed in the seminal work of Meier and Zierler22 and then extended to the general case of extravascular tracers (as in the case of PET) by O'Sullivan et al,3 an equation analog to (1) can be written: Z t C T ðtÞ ¼ C P IRFðtÞ ¼ K Rðt - sÞC P ðsÞds ð5Þ 0
Equation (5) provides a probabilistic framework for representation of PET tracer data in terms of a convolution between the metabolite-corrected input function and a tissue residue function R(t). K is a proportionality constant, interpretable as an overall flow. The residue is a scaled survival Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1370 function associated with tracer residence in the tissue. Physically, the residue function describes the fraction of indicator that remains in the tissue in response to an idealized bolus (Dirac delta function) input concentration at time t = 0. Initially, the residue must be unity (R(0) = 1) and from there it decreases (at least does not increase) as t increases. Assuming the residue is differentiable, it can be represented in the form: Z t hðτÞdτ ð6Þ RðtÞ ¼ 1 0
where h(t) is a probability density function. Meier and Zierler22 developed an interpretation of h(t) in terms of the residence time distribution associated with an ensemble of capillary transport paths within the region. Indicator molecules are assumed to randomly sample such paths. Once a path is selected, the travel time of the molecule in the tissue is determined. Either due to flow, volume variations or binding, the transit time of the indicator in each transport path is considered to be distinct. From the distribution of transport paths there arises a distribution of transit times that would be experienced by a randomly chosen indicator molecule. The probability that an indicator molecule ends up on a path with transit time between t and t+Δt is proportional to h(t)Δt as Δt gets small. Hence if a collection of indicator molecules is introduced in trace amounts, so they do not interact, (6) gives the fraction of molecules that experiences transit times greater than t.
Outcome Measures Expressing (1) as in (5), and considering nonparametric inference for the residue, gives rise to a deconvolution problem, providing a novel approach to the analysis of dynamic PET data that is not reliant on specific compartmental modeling assumptions.3 However, once the IRF has been nonparametrically estimated via deconvolution, it remains to establish which are the outcomes that can be actually computed from the IRF, and how they compare with respect to traditional outcomes used in PET. The first outcome measure usually computed in PET for a reversible tracer is VT, whose formulation in terms of the 2TCM K1 to k4 values is:2 Z 1 K1 k3 VT ¼ 1þ f ðt; k 2 ; k 3 ; k 4 Þdt ð7Þ ¼ K1 k2 k4 0 VT equals the integral up to infinity of the corresponding IRF, which is only made possible by the parametric formulation of the IRF. When the IRF is computed nonparametrically, however, the IRF is estimated only up to the end time of the PET scan. A different definition is therefore necessary, as proposed by Jiang et al,4 in which a surrogate ‘volume of distribution’ for the tracer is calculated by integrating the IRF only up to the end time Tend of the scan: Z T end V T - end ¼ K RðπÞdπ ≠V T ð8Þ 0
To obtain a volume that is comparable to the approach used with CMs, O’Sullivan et al3 proposed an alternative formulation for VT based on the idea that, if the residence density is partitioned into quantiles, the volumes occupied by fast, intermediate, and slow moving tracer label can be separately assessed.3 Their formulation uses the scaled first quantile of the residence density: V T - NP1 ¼
KR - 1 ð1=4Þ log ð4Þ
−1
ð9Þ
where R (1/4) indicates the first quantile of the residence, which is the time at which the IRF reaches a fourth of its initial value, and the scaling log(4) ensures that if the residence density is exponential (i.e., as in the 1TCM case), then R − 1(1/4)/log(4) is equal to the mean of the residence distribution. The scaling factor in Equation (9) differs from the one reported in the paper by O’Sullivan et al.3 We closely checked our derivation of the scaling factor under the assumption that the residence density is exponential (h(t) = α × exp(− αt)). We confirmed our derivation using the reported clinical and simulated data. A closed-form expression for the scaling factor when the residence density is biexponential or more cannot be derived. While the use of the first quantile of the residence was suggested for the irreversible tracer [18F]FDG,3 we propose here to further generalize and optimize this nonparametric formulation of VT for reversible tracers, by considering not only the first quantile of the residence density (i.e., R −1(1/4)), but also a large number of quantiles (i.e., R − 1(n), for n = N1, …, NQ), Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
and then averaging over all of them: V T - NP2 ¼
K Xn¼NQ R - 1 ðnÞ n¼N1 log ð1=nÞ N
ð10Þ
This is a natural extension of the idea of tying the nonparametric estimate of the IRF with a particular counterpart; rather than relying only on a single estimated quantile with its inherent uncertainty, this approach will tend to ‘sharpen’ the estimator by considering many estimated quantiles with the expectation that the uncertainties will cancel each other out to some extent. Usually with CMs, once VT is computed according to (7), the binding potential BPP in any ROI/voxel can be indirectly estimated as BPP = VT − VTref, where VTref is the tracer VT in an RR devoid of the target of interest.1 The same formulation of BPP can be adopted for the VT calculated according to Equations (8, 9, and 10), to obtain BPP-end, BPP-NP1, and BPP-NP2, respectively.
Singular Value Decomposition Potentially, any of many nonparametric deconvolution approaches can be used to estimate the IRF starting from known metabolite-corrected input function and tracer TACs in the brain tissue. O’Sullivan et al3 propose a technique based on regularized cubic B-spline approximation of the residence time distribution, but do not claim any optimality of the B-spline estimators. Jiang et al4 use functional principal component analysis. Here, we propose the use of SVD. Specifically, if the tissue concentration CT(t) of the tracer in response to a metabolite-corrected input function CP(t) is given by Equation (5), then assuming that the metabolite-corrected input function and the tracer TACs in the brain tissue are measured at equally spaced time points t1, t2, …, tN, with Δt as the space between time points of measurement, Equation (5) can be discretized, assuming that over short time intervals Δt, the residue function and metabolite-corrected input function are constant in time: Z tj Xj C P ðπÞRðt - πÞdπ ffi KΔt C ðt ÞRðtj - t i Þ ð11Þ C T ðtj Þ ¼ K i¼0 P j 0
or
0
C P ðt1 Þ B C P ðt2 Þ B KΔt@ … C P ðtN Þ
0 C P ðt1 Þ … C P ðtN - 1 Þ
10 1 0 1 Rðt1 Þ C T ðt1 Þ … 0 B Rðt2 Þ C B C T ðt2 Þ C … 0 C CB C¼B C … … A@ … A @ … A RðtN Þ C T ðtN Þ … C P ðt1 Þ
ð12Þ
The requirement for equidistantly sampled data derives from the Shannon’s sampling theorem,23 which implies that any band-limited signal can be numerically reconstructed from its regularly spaced samples. The solution of the nonuniform sampling problem involves much more complex theoretical and computational issues (see Discussion).24 In matrix form, (12) is expressed as: Ab ¼ c
ð13Þ
where b contains the elements of R(ti), i = 1, 2,…, N, and c are the measured tracer TACs in the brain tissue CT(ti), i = 1, 2,…, N. Equation (13) can be solved iteratively for the elements of b. This approach is, however, extremely sensitive to noise, potentially causing the estimated R(t) to oscillate. A technique to solve (13) is SVD, which constructs matrices V, W, and UT, so that the inverse of A in (13) can be written: A - 1 ¼ VWUT
ð14Þ T
where W is a diagonal matrix and V and U are orthogonal matrices. Given this inverse matrix, b, and consequently R(t), is calculated as b ¼ A - 1 c ¼ VWUT c
ð15Þ
The main force of the SVD is that the diagonal elements in W are zero or close to zero corresponding to linear equations in (12) that are close to being linear combinations of each other. This fact thus allows one to identify elements in the matrix A that causes the solution b to oscillate or otherwise be meaningless in a biomedical modeling context. By eliminating (setting equal to zero) diagonal elements below a certain threshold in W, one can consequently minimize these effects before calculating b.25 Here, we propose an automatic data-driven algorithm for selecting the threshold for SVD for each ROI that is based on cross-validation. Specifically, for each TAC in each considered ROI, the threshold is automatically selected by: © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1371 1. For i = 1 to i = N (with N the total number of original sampling time points) (1) Removing the TAC point corresponding to the sampling time point ti; (2) For threshold T = 2% to T = 40% (with step = 1%), applying SVD with threshold value T to the TAC containing the remaining N − 1 data points; (3) Storing the threshold T that gives the best prediction of the TAC data at the time point ti that was left out (i.e., calculated as the minimum of the square difference between original TAC value at the excluded time point ti and the reconvolved TAC obtained with SVD at the same ti); 2. Selecting the average value (across excluded time points ti) of the stored thresholds. The suggested threshold range is a generous range that has proven optimal before for physiologic signals.9,10
Metrics of Comparison for Simulated Data For all simulated cases and noise conditions, we ran SVD with the automatic threshold selection, and calculated VT-NP1, VT-NP2, and VT-end from the deconvolved IRF estimates. BPP-NP1, BPP-NP2, and BPP-end were derived by subtracting VT-NP1, VT-NP2, and VT-end of the RR. We compared the distribution of outcome measures estimates obtained with the true simulated values, and calculated variance (across the 1,000 estimates of the outcome measure in each ROI, the median variance across ROIs is reported) and bias (as distance between median of the 1,000 estimates of the outcome measure in each ROI and the corresponding true value, the median bias across ROIs is then considered) of such outcome measures.
Application to Test–Retest Data For all tracers, subjects, and ROIs, we calculated VT and BPP = VT − VTref using CMs (2TCM for CUMI,14 PE2I,16 and WAY;17 1TCM for DASB15); the RR was the gray matter (for CUMI14 and DASB15), the white matter (for WAY17), or the totality (for PE2I16) of the cerebellum. We also ran SVD with the automatic threshold selection, and calculated VT-NP1, VT-NP2, and VT-end from the deconvolved IRF. BPP-NP1, BPP-NP2, and BPP-end were derived by subtracting VT-NP1, VT-NP2, and VT-end, respectively, of the RR. We first compared VT-NP1, VT-NP2, and VT-end versus VT values using correlation (r), slope and intercept of the regression line, with VT values being the independent variables. BPP-NP1, BPP-NP2, and BPP-end were similarly compared with BPP values. We then compared the test–retest performance of VT, VT-NP1, VT-NP2, VT-end, and BPP, BPP-NP1, BPP-NP2, BPP-end, and using five metrics:15 (1) Percent difference (PD), calculated as the absolute difference between the test and retest values divided by their average (the lowest the PD, the highest the reproducibility); (2) Within-subject of squares (WSMSS): sum P P mean WSMSS ¼ n1 ni¼1 2j¼1 X ij - X i , where Xi1 represents the ‘test’ value of an outcome measure for subject i, Xi2 the ‘retest’ value for the same subject, n is the number of pairs of subjects, and X i ¼ ðX i1 þ X i2 Þ=2. This represents the variance in the estimation of a subject-specific outcome measure (the lowest the WSMSS, the highest the reproducibility); BSMSS - WSMSS (3) Interclass correlation coefficient (ICC): ICC ¼ BSMSSþðk - 1ÞWSMSS, where BSMSS is the between-subject mean sum of squares
P 1 Pk and k is the number of repeated X X BSMSS ¼ n1 2n ij i i¼1 k j¼1 observations (k = 2 in this study; − 1 corresponds to no reliability, 1 to maximum reliability); (4) Variance of an outcome measure across subjects (sample variance, across subjects, of all observations for each ROI); and (5) Identifiability based on bootstrap resampling of residuals; we computed 100 bootstrap samples26 of the TAC data for each subject’s ROI, estimated the outcome measures for each of these samples, and measured the variability of these estimates using the robust median absolute deviation criterion: median [|Zj – median(Z1, Z2, y , Z100)|, j = 1,2, y ,100]. © 2015 ISCBFM
RESULTS Test–Retest Data Comparison with CMs estimates. Figures 1 and 2 show for all tracers, subjects, and ROIs, the scatter plots of VT-NP1, VT-NP2, and VT-end, and BPP-NP1, BPP-NP2, and BPP-end, respectively, obtained using SVD versus the VT and BPP values obtained using CMs. The identity and regression lines, and the numerical results of the regression analysis, are reported in each plot. Correlation is high (40.90) in all cases, with exception of BPP-NP1 for PE2I (r = 0.740) and BPP-end for CUMI (r = 0.875). However, the slope values differ across outcome measures and tracers. As expected from its formulation in comparison with VT, VT-end underestimation increases as the value of the estimated outcome increases. VT-end shows the lowest slope values across the different outcomes for all tracers: the minimum value is 0.466 (PE2I) and the maximum is 0.619 (DASB). VT-NP1 and VT-NP2 show similar results with each other, and the slope values are closer to unity than with VT-end (VT-NP1: minimum value is 0.511, PE2I, maximum is 1.019, WAY; VT-NP2: minimum value is 0.827, PE2I, maximum is 0.944, DASB). As expected by definition, VT-NP2 shows higher correlation than VT-NP1 for all tracers, while the slope is closer to unity than with VT-NP1 for DASB and PE2I. For all tracers, slope values for BPP decrease with respect to the comparison with VT values. In general, BPP-NP2 shows the best correlation values for all tracers (0.922 for CUMI, 0.986 for DASB, 0.973 for PE2I, and 0.942 for WAY), and the best slope values for DASB (0.928) and PE2I (0.794). BPP-NP1 shows the best slope value for CUMI (0.858) and WAY (1.044). Test–retest metrics. Figures 3 and 4 show the boxplots of the five considered metrics for all tracers and formulations of the outcome measure (VT, VT-NP1, VT-NP2, VT-end and BPP, BPP-NP1, BPP-NP2, BPP-end, respectively). In each box, the central mark is the median value of the metric (across subjects/pairs and ROIs), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. A complete list of the numerical results corresponding to Figures 3 and 4 is reported in the Supplementary Material. Table 1 summarizes which outcome measure provides the best performance. For the clinical data, the outcome measures with the best correlation and slope (with the corresponding outcome measures estimated with CMs) are reported separately from the test–retest metrics; the best performance in the test–retest metrics is reported both considering and not considering VT-end (or BPP-end). VT-end shows the lowest WSMSS and variance, and the best identifiability, for all tracers. VT-end also shows the lowest PD for CUMI, DASB, and PE2I, and the highest ICC for CUMI and DASB. Correspondingly, BPP-end shows the lowest WSMSS and variance for all tracers. BPP-end also shows the highest ICC for all tracers, besides WAY, the lowest PD for DASB and PE2I, and the best identifiability for all tracers, with exclusion of DASB. VT and BPP show the lowest PD and the highest ICC for WAY. BPP also shows the best identifiability for DASB. VT-NP1 and BPP-NP1 show the worst identifiability for all tracers. BPP-NP1 also shows the lowest ICC for all tracers, besides WAY. VT-NP2 shows the highest ICC for PE2I (0.80), while BPP-NP2 shows the lowest PD for CUMI. Simulated Data Figures 5 and 6 show the distribution of outcome measures obtained using SVD in the 1TCM and 2TCM case, respectively (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2), with simulated noise level corresponding to α = 1. Each region is Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1372
Figure 1. Scatter plots of VT-NP1, VT-NP2, and VT-end obtained by singular value decomposition (SVD) versus the VT values obtained using compartment models (CMs) for all considered tracers, subjects, and regions of interest (ROIs). The identity (solid) and regression (dotted) lines are reported in each plot, together with the numerical results of the regression analysis. VT, total volume of distribution.
reported separately (the cerebellum gray matter is the RR for calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation of the outcome measure, while the black solid vertical line represents the true simulated value for that outcome measure and region. In the 1TCM case, VT-end always underestimates VT the most across the three formulations, except for the cerebellum gray matter, although with much less variability. In the 2TCM case, VT-end performs the best for the three regions with the lowest VT (cerebellum gray matter, occipital lobe, and cingulate), but underestimates the most the region with the highest VT (hippocampus). Correspondingly, BPP-end always underestimates BPP the most in both cases. In the 1TCM case, VT-NP1 and VT-NP2 show similar median values for VT, with VT-NP2 performing better than VT-NP1 for the three Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
regions with highest simulated VT (dorsal caudate, amydgala, ventral striatum). In the 2TCM case, it is VT-NP1 that slightly outperforms VT-NP2 in all regions in terms of bias, but not variance of the estimates distribution. Corresponding BPP-NP1 and BPP-NP2 values are similar in the 1TCM case, while BPP-NP1 tends to overestimate the BPP in all the regions in the 2TCM case. BPP-NP2 shows the smallest bias across regions and less dispersed distributions than the ones of BPP-NP1 in the 2TCM case. As reported in Table 1, with the noise present at the ROI level (α = 1), VT-end and BPP-end show the best variance in both simulated cases: 0.039 and 0.084 (1TCM); 0.043 and 0.090 (2TCM). VT-end shows the smallest bias (−0.298) for 2TCM, while VT-NP2 shows the smallest bias for 1TCM (−0.426), and BPP-NP2 shows the smallest bias for both 1TCM (0.129) and 2TCM (0.351). The results relative to the noise analysis are reported in the Supplementary Material. In both 1TCM and 2TCM cases, the bias of © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1373
Figure 2. Scatter plots of BPP-NP1, BPP-NP2, and BPP-end obtained by singular value decomposition (SVD) versus the BPP values obtained using compartment models (CMs) for all considered tracers, subjects, and regions of interest (ROIs). The identity (solid) and regression (dotted) lines are reported in each plot, together with the numerical results of the regression analysis. BPp, binding potential.
the VT-end and BPP-end estimates is almost constant across the five noise conditions. The corresponding variance, although it increases with the noise in the data, remains significantly lower than the corresponding variance of the other two nonparametric formulations for the same outcome measure. In both 1TCM and 2TCM cases, the pattern with which the bias of VT-NP1 and BPP-NP1, and VT-NP2 and BPP-NP2, changes with noise is similar, with both formulations that tend to overestimate the true value at elevated noise levels. In the 1TCM case, the bias of VT-NP2 and BPP-NP2 increases the most with the noise level. The corresponding variance values increase with noise, with VT-NP1 and BPP-NP1 providing the worst performance in the 2TCM case. In the 1TCM case, the variance of VT-NP1 and VT-NP2, and BPP-NP1 and BPP-NP2, is almost identical. © 2015 ISCBFM
DISCUSSION We applied nonparametric deconvolution to estimate the IRF of PET reversible tracers in several regions, using both simulated and clinical test–retest data sets, computed several functionals of the IRF to estimate tracer VT and BPs, and extensively compared these functionals with traditional outcome measures derived using CMs. We showed that nonparametric deconvolution, which does not require any model assumptions, and the functionals that can be computed from the recovered IRF can provide estimates of VT and BPP that are not only very close to, and highly correlated with, the estimated obtained with traditional CMs, but can also show better test–retest performance than outcome measures derived from CMs. Therefore, nonparametric deconvolution and the proposed functionals offer a promising alternative to CMs for the Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1374
Figure 3. Boxplots of the five considered test–retest metrics for all considered tracers and the four different estimations of the outcome measure (VT, VT-NP1, VT-NP2, VT-end). In each box, the central mark is the median value of the metric (across subjects/pairs and regions), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. ICC, interclass correlation coefficient; PD, percent difference; VT, total volume of distribution; WSMSS, within-subject mean sum of squares.
quantification of PET dynamic data, especially when CMs assumptions do not hold or when CMs provide poor fits and/or unphysiologic estimates of the outcome measures for a specific PET tracer. The optimized formulation of VT we proposed in Equation (10), and the corresponding BPP-NP2 that can be derived, provide the best performance in terms of bias with respect to, and correlation with, outcome measures with traditional CMs in the majority of the simulated and clinical considered data sets (Table 1). VT-end as in Equation (8), and the corresponding BPP-end that can be derived, provide the best test–retest performance across metrics and, overall, considered tracers; outcome measures derived with traditional CMs show superior test–retest performance than VT-end and BPP-end only for some of the metrics in the case of DASB Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
(identifiability of BPP), PE2I (ICC of VT), and WAY (PD and ICC of VT and BPP). Therefore, VT-end and BPP-end are stable (with respect to noise in the TACs) and reliable outcome measures that could be used, for example, to discriminate different clinical populations from each other and from healthy controls. However, as one would expect from their definitions, results on simulated and clinical data sets suggest that VT-end and BPP-end often provide the greatest underestimation of VT and BPP (with the exception of PE2I BPP, for which BPP-NP1 shows the greatest underestimation). This is perhaps less of an issue if the purpose of the study is to examine differences between groups (or across values of some continuous variable) but if the primary goal of the analysis is estimation of standard outcome measures, the bias of these estimates must be considered. In such a situation, VT-NP1 and/or VT-NP2 (and the © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1375
Figure 4. Boxplots of the five considered test–retest metrics for all considered tracers and the four different estimations of the outcome measure (BPP, BPP-NP1, BPP-NP2, BPP-end). In each box, the central mark is the median value of the metric (across subjects/pairs and regions), the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points of the metric not considered outliers, and outliers are plotted individually. BP, binding potential; ICC, interclass correlation coefficient; PD, percent difference; WSMSS, within-subject mean sum of squares.
corresponding BPP-NP1 and BPP-NP2) could be considered, as they show better test–retest performance than CMs outcomes according to the majority of the considered metrics for all tracers, with the exception of WAY. An intrinsic limitation of the nonparametric approach is that the IRF can be estimated only up to the end of the PET scan, and no inference up to infinity can therefore be made without some strong assumptions. For example, it may be the case that, for tracers that slowly reach equilibrium, the nonparametric estimation of VT and BPP according to Equations (9 and 10) is biased by the fact that, by the end of the scan, the IRF has not yet reached a certain fraction of its initial value, and therefore the corresponding kth quantile of the residence cannot be determined. Estimation of VT-end and BPP-end instead might be preferred in these cases, with proper interpretation of the results. © 2015 ISCBFM
We have implemented here nonparametric deconvolution using SVD9,10 with an automatic data-driven selection of the threshold for each ROI. This is a fast, easy to implement, and widely adopted approach to deconvolution that easily lends itself to implementation at the voxel level. However, we claim no optimality in either the choice of the deconvolution approach or the selection of the thresholds. It is known that SVD performance can be sensitive to the level of the noise present in the data.27,28 Our simulation results at different noise level (see Supplementary material) suggest that, at the voxel level, or for particularly noisy tracers, either VT-end and BPP-end should be preferred as outcome measures because they are less affected by noise, or more robust approaches to deconvolution, such as nonlinear stochastic regularization28 or Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1376 Table 1.
Summary of which outcome measure provides the best performance in each of the simulated and clinical data sets Simulated data
1TC 2TC Clinical data CUMI DASB PE2I WAY
VT Best variance
Best bias
BPp Best variance
Best bias
VT-end VT-end
VT-NP2 VT-end
BPp-end BPp-end
BPp-NP2 BPp-NP2
Best correlation
Best slope
Best correlation
Best slope
VT-NP2 VT-NP2 VT-NP2 VT-NP2
VT-NP1 VT-NP2 VT-NP2 VT-NP2
BPp-NP2 BPp-NP2 BPp-NP2 BPp-NP2
BPp-NP1 BPp-NP2 BPp-NP2 BPp-NP2
Test–retest metrics
VT (considering VT-end) CUMI DASB PE2I WAY VT (not considering VT-end) CUMI DASB PE2I WAY BPp (considering BPp-end) CUMI DASB PE2I WAY BPp (not considering BPp-end) CUMI DASB PE2I WAY
PD
WSMSS
ICC
Variance
Identifiability
VT-end VT-end VT-end VT
VT-end VT-end VT-end VT-end
VT-end VT-end VT/VT-NP2 VT
VT-end VT-end VT-end VT-end
VT-end VT-end VT-end VT-end
VT-NP2 VT VT-NP1 VT
VT-NP2 VT-NP2 VT-NP1 VT
VT-NP2 VT/VT-NP2 VT/VT-NP2 VT
VT-NP2 VT-NP1 VT-NP1 VT
VT-NP2 VT VT-NP2 VT
BPp-NP2 BPp-end BPp-end BPp
BPp-end BPp-end BPp-end BPp-end
BPp-end BPp-end BPp-end BPp
BPp-end BPp-end BPp-end BPp-end
BPp-end BPp BPp-end BPp-end
BPp-NP2 BPp BPp BPp
BPp-NP2 BPp-NP2 BPp-NP1 BPp
BPp-NP2 BPp-NP2 BPp-NP2 BPp
BPp-NP2 BPp BPp-NP1 BPp
BPp BPp BPp-NP2 BPp
Abbreviations: BPp, binding potential; CM, compartment model; ICC, interclass correlation coefficient; PD, percent difference; 1TC, one-tissue compartment; 2TC, two-tissue compartment; VT, total volume of distribution; WSMSS, within-subject mean sum of squares. For the simulated data, the outcome measures with the best variance and best bias (with respect to the true value of the outcome measure) are reported. For the clinical data, the outcome measures with the best correlation and slope (with the corresponding outcome measures estimated with CM) are reported separately from the test–retest metrics; the best performance in the test–retest metrics is reported considering and not considering VT-end (or BPP-end).
kernel-based linear system identification,29,30 should be considered for reconstructing the IRF. The cross-validation algorithm we propose here for thresholding the singular values in the decomposition to constrain the propagation of noise from the data to the IRF estimates is data driven and has proven to be flexible enough to provide reasonably regularized solutions for the considered four tracers and simulated data. However, other strategies exist,7,31–33 which may provide improved solutions. The optimization of the thresholding algorithm, which likely depends on the tracer at hand, is however beyond the scope of this article. Singular value decomposition, as many other nonparametric deconvolution approaches, requires the input data to be equidistantly sampled.23 The PET data involve radioactive decay, and are usually acquired using increasing frame lengths to ensure a comparable number of counts and thus statistical quality per data point. This implies that the PET data need to be interpolated on an equally-spaced grid of temporal values before applying SVD. The possible effect of the interpolation on the estimation of VT and BPP should be further investigated, especially when using SVD at the voxel level where, due to the large number of voxels to quantify, adopting a very fine interpolation grid can Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
lead to out-of-memory issues during the quantification process. While Shannon’s sampling theorem23 implies that any band-limited signal can be reconstructed from its regularly spaced samples, the solution of the nonuniform sampling problem poses more difficulties.24 There exist a vast literature on nonuniform sampling theory34,35 and methods claiming to efficiently reconstruct a function from its nonuniform samples,36–38 including nonuniform SVD.39 The investigation of the performance on PET data of nonuniform deconvolution approaches is however beyond the scope of this article. It is important to state that comparing the performance of the proposed nonparametric outcome measures with the corresponding CMs estimates is not generally useful when the CM is not valid for the tracer in question. Here, we have considered four tracers for which we have already assessed the validity of the CMs,14–17 so the proposed comparison is meaningful, although it may be questioned in certain individual regions. Certainly, the possibility of quantifying the binding of PET tracers model-free and nonparametrically is particularly attractive for those tracers for which the assumptions of CMs do not hold or CMs do not provide a good description of the data. In these cases, however, the CMs estimates of the binding cannot be used as a © 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1377
Figure 5. Distribution of outcome measures estimated using singular value decomposition (SVD) and the three nonparametric formulations (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2) across the 1,000 simulated one-tissue compartment model (1TCM) cases in each region; the six simulated regions are reported separately (the cerebellum gray matter is used as the reference region (RR) for the calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation, while the black solid vertical line represents the true simulated value for that outcome measure and region. The data underlying the simulations were derived from clinical [11C]DASB scans modeled with 1TCM. BPp, binding potential; VT, total volume of distribution.
reference to evaluate the performance of the nonparametric approach, and validation with in vitro results would be recommended.
© 2015 ISCBFM
Future directions of investigation also include the extension of the presented nonparametric approach to quantification in the absence of a metabolite-corrected input function.
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
Model-free nonparametric PET quantification F Zanderigo et al
1378
Figure 6. Distribution of outcome measures estimated using singular value decomposition (SVD) and the three nonparametric formulations (top panel: VT-end, VT-NP1, VT-NP2; bottom panel: BPP-end, BPP-NP1, BPP-NP2) across the 1,000 simulated two-tissue compartment model (2TCM) cases in each region; the five simulated regions are reported separately (the cerebellum gray matter is used as the reference region (RR) for the calculation of BPP-end, BPP-NP1, BPP-NP2; therefore, it is not reported for that outcome measure); the dotted vertical lines represent the median of the distribution for each formulation, while the black solid vertical line represents the true simulated value for that outcome measure and region. The data underlying the simulations were derived from clinical [11C]CUMI scans modeled with 2TCM. BPp, binding potential; VT, total volume of distribution.
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
© 2015 ISCBFM
Model-free nonparametric PET quantification F Zanderigo et al
1379 CONCLUSIONS Nonparametric deconvolution does not require any model assumptions and allows for estimates of VT and BPP that are very close to the estimates obtained with traditional CMs, and can also show superior test–retest performance. Nonparametric deconvolution and the proposed functionals therefore offer a viable alternative to CMs for the quantification of PET dynamic data, especially when CMs assumptions do not hold or CMs provide poor fits and/or unphysiologic estimates of the outcome measures for a specific PET tracer. AUTHOR CONTRIBUTIONS Dr Zanderigo was responsible for ideation of the methodologies, drafting of the manuscript, creation of the software, analysis and interpretation of the data. Dr Parsey contributed to conception and revision of the manuscript, and interpretation of the data. Dr Ogden contributed to ideation of the methodologies, interpretation of the data, critical revision of the manuscript, and final approval.
DISCLOSURE/CONFLICT OF INTEREST The authors declare no conflict of interest.
REFERENCES 1 Innis RB, Cunningham VJ, Delforge J, Fujita M, Gjedde A, Gunn RN et al. Consensus nomenclature for in vivo imaging of reversibly binding radioligands. J Cereb Blood Flow Metab2007; 27: 1533–1539. 2 Slifstein M, Laruelle M. Models and methods for derivation of in vivo neuroreceptor parameters with PET and SPECT reversible radiotracers. Nucl Med Biol2001; 28: 595–608. 3 O'Sullivan F, Muzi M, Spence AM, Mankoff DM, O'Sullivan JN, Fitzgerald N et al. Nonparametric residue analysis of dynamic PET data with application to cerebral FDG studies in normals. J Am Stat Assoc 2009; 104: 556–571. 4 Jiang C-R, Aston JAD, Wang J-L. A functional approach to deconvolve dynamic neuroimaging data 2014; arXiv.org/stat/arXiv:1411.2051. 5 Cunningham VJ, Jones T. Spectral analysis of dynamic PET studies. J Cereb Blood Flow Metab 1993; 13: 15–23. 6 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: 1425–1439. 7 Turkheimer FE, Hinz R, Gunn RN, Aston JA, Gunn SR, Cunningham VJ. Rankshaping regularization of exponential spectral analysis for application to functional parametric mapping. Phys Med Biol 2003; 48: 3819–3841. 8 Schmidt K. Which linear compartmental systems can be analyzed by spectral analysis of PET output data summed over all compartments? J Cereb Blood Flow Metab 1999; 19: 560–569. 9 Ostergaard L, Weisskoff RM, Chesler DA, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis. Magn Reson Med 1996; 36: 715–725. 10 Ostergaard L, Sorensen AG, Kwong KK, Weisskoff RM, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part II: Experimental comparison and preliminary results. Magn Reson Med 1996; 36: 726–736. 11 Buchert R, Wilke F, van den Hoff J, Mester J. Improved statistical power of the multilinear reference tissue approach to the quantification of neuroreceptor ligand binding by regularization. J Cereb Blood Flow Metab 2003; 23: 612–620. 12 Constantinescu CC, Bouman C, Morris ED. Nonparametric extraction of transient changes in neurotransmitter concentration from dynamic PET data. IEEE Trans Med Imaging 2007; 26: 359–373. 13 Grecchi E, Doyle OM, Bertoldo A, Pavese N, Turkheimer FE. Brain shaving: adaptive detection for brain PET data. Phys Med Biol 2014; 59: 2517–2534. 14 Milak MS, DeLorenzo C, Zanderigo F, Prabhakaran J, Kumar JS, Majo VJ et al. In vivo quantification of human serotonin 1A receptor using 11C-CUMI-101, an agonist PET radiotracer. J Nucl Med 2010; 51: 1892–1900.
15 Ogden RT, Ojha A, Erlandsson K, Oquendo MA, Mann JJ, Parsey RV. In vivo quantification of serotonin transporters using [(11)C]DASB and positron emission tomography in humans: modeling considerations. J Cereb Blood Flow Metab 2007; 27: 205–217. 16 DeLorenzo C, Kumar JS, Zanderigo F, Mann JJ, Parsey RV. Modeling considerations for in vivo quantification of the dopamine transporter using [(11)C]PE2I and positron emission tomography. J Cereb Blood Flow Metab 2009; 29: 1332–1345. 17 Parsey RV, Slifstein M, Hwang DR, Abi-Dargham A, Simpson N, Mawlawi O et al. Validation and reproducibility of measurement of 5-HT1A receptor parameters with [carbonyl-11C]WAY-100635 in humans: comparison of arterial and reference tisssue input functions. J Cereb Blood Flow Metab 2000; 20: 1111–1133. 18 Jenkinson M, Smith S. A global optimisation method for robust affine registration of brain images. Med Image Anal 2001; 5: 143–156. 19 Smith SM. Fast robust automated brain extraction. Hum Brain Mapp 2002; 17: 143–155. 20 Parsey RV, Ogden RT, Mann JJ. Determination of volume of distribution using likelihood estimation in graphical analysis: elimination of estimation bias. J Cereb Blood Flow Metab 2003; 23: 1471–1478. 21 Parsey RV, Ojha A, Ogden RT, Erlandsson K, Kumar D, Landgrebe M et al. Metabolite considerations in the in vivo quantification of serotonin transporters using 11C-DASB and PET in humans. J Nucl Med 2006; 47: 1796–1802. 22 Meier P, Zierler KL. On the theory of the indicator-dilution method for measurement of blood flow and volume. J Appl Physiol 1954; 6: 731–744. 23 Shannon C. A mathematical theory of communication. Bell Syst Tech J 1948; 27: 379–623. 24 Strohmer T. Numerical analysis of the non-uniform sampling problem. J Comput Appl Math 2000; 122: 297–316. 25 Press WH, Teukolsky SA, Vetterling WT, Flannery BT. Numerical Recipes in C. 2nd edition Cambridge University Press: Oxford, 1992. 26 Ogden RT, Tarpey T. Estimation in regression models with externally estimated parameters. Biostatistics 2006; 7: 115–129. 27 Knutsson L, Stahlberg F, Wirestam R. Aspects on the accuracy of cerebral perfusion parameters obtained by dynamic susceptibility contrast MRI: a simulation study. Magn Reson Imaging 2004; 22: 789–798. 28 Zanderigo F, Bertoldo A, Pillonetto G, Cobelli Ast C. Nonlinear stochastic regularization to characterize tissue residue function in bolus-tracking MRI: assessment and comparison with SVD, block-circulant SVD, and Tikhonov. IEEE Trans Biomed Eng 2009; 56: 1287–1297. 29 Pillonetto G. DNG. A new kernel-based approach for linear system identification. Automatica 2010; 46: 81–93. 30 Castellaro M, Peruzzo D, Mehndiratta A, Pillonetto G, Petersen ET, Golay X et al. Estimation of arterial arrival time and cerebral blood flow from QUASAR arterial spin labeling using stable spline. Magn Reson Med 2014; doi: 10.1002/mrm.25525 (e-pub ahead of print). 31 Murase K, Yamazaki Y, Miyazaki S. Deconvolution analysis of dynamic contrastenhanced data based on singular value decomposition optimized by generalized cross validation. Magn Reson Med Sci2004; 3: 165–175. 32 Sourbron S, Luypaert R, Van Schuerbeek P, Dujardin M, Stadnik T. Choice of the regularization parameter for perfusion quantification with MRI. Phys Med Biol2004; 49: 3307–3324. 33 Liu HL, Pu Y, Liu Y, Nickerson L, Andrews T, Fox PT et al. Cerebral blood flow measurement by dynamic contrast MRI using singular value decomposition with an adaptive threshold. Magn Reson Med 1999; 42: 167–172. 34 Landau HJ. Necessary density conditions for sampling and interpolation of certain entire functions. Acta Math 1967; 117: 37–52. 35 Feichtinger HG, Gröchenig KH. Theory and practice of irregular sampling. Benedetto J, Frazier M, (Eds.), Wavelets: Mathematics and Applications. CRC Press: Boca Raton, FL, 1994; 305–363. 36 Yao K TJ. On some stability and interpolatory properties of nonuniform sampling expansions. IEEE Trans Circuit Theory 1967; 14: 404–408. 37 Feichtinger HG GK, Strohmer T. Efficient numerical methods in non-uniform sampling theory. Numer Math 1995; 69: 423–440. 38 De Nicolao G, Liberati D, Sartorio A. Deconvolution of infrequently sampled data for the estimation of growth hormone secretion. IEEE Trans Biomed Eng 1995; 42: 678–687. 39 Wingham D. The reconstruction of a band-limited function and its Fourier transform from a finite number of samples at arbitrary locations by singular value decomposition. IEEE Trans Circuit Theory 1967; 40: 559–570.
Supplementary Information accompanies the paper on the Journal of Cerebral Blood Flow & Metabolism website (http://www.nature. com/jcbfm)
© 2015 ISCBFM
Journal of Cerebral Blood Flow & Metabolism (2015), 1368 – 1379
