A hybrid deconvolution approach for estimation of in vivo non-displaceable binding for brain PET targets without a reference region

Background and aim Estimation of a PET tracer’s non-displaceable distribution volume (VND) is required for quantification of specific binding to its target of interest. VND is generally assumed to be comparable brain-wide and is determined either from a reference region devoid of the target, often not available for many tracers and targets, or by imaging each subject before and after blocking the target with another molecule that has high affinity for the target, which is cumbersome and involves additional radiation exposure. Here we propose, and validate for the tracers [11C]DASB and [11C]CUMI-101, a new data-driven hybrid deconvolution approach (HYDECA) that determines VND at the individual level without requiring either a reference region or a blocking study. Methods HYDECA requires the tracer metabolite-corrected concentration curve in blood plasma and uses a singular value decomposition to estimate the impulse response function across several brain regions from measured time activity curves. HYDECA decomposes each region’s impulse response function into the sum of a parametric non-displaceable component, which is a function of VND, assumed common across regions, and a nonparametric specific component. These two components differentially contribute to each impulse response function. Different regions show different contributions of the two components, and HYDECA examines data across regions to find a suitable common VND. HYDECA implementation requires determination of two tuning parameters, and we propose two strategies for objectively selecting these parameters for a given tracer: using data from blocking studies, and realistic simulations of the tracer. Using available test-retest data, we compare HYDECA estimates of VND and binding potentials to those obtained based on VND estimated using a purported reference region. Results For [11C]DASB and [11C]CUMI-101, we find that regardless of the strategy used to optimize the tuning parameters, HYDECA provides considerably less biased estimates of VND than those obtained, as is commonly done, using a non-ideal reference region. HYDECA test-retest reproducibility is comparable to that obtained using a VND determined from a non-ideal reference region, when considering the binding potentials BPP and BPND. Conclusions HYDECA can provide subject-specific estimates of VND without requiring a blocking study for tracers and targets for which a valid reference region does not exist.


Introduction
Positron Emission Tomography (PET) in the brain involves administration of a tracer dose of a radioactively labeled molecule (i.e., tracer) that binds to a specific target [1].The tracer signal in the tissue combines signal from tracer "specifically" bound to the target and tracer "nonspecifically" bound to other macromolecules or free in tissue water.Estimation of tracer nondisplaceable uptake allows quantification of the specific binding potential between tracer and target [2,3].The tracer non-displaceable distribution volume (V ND ), corresponding to "nonspecifically" bound and free tracer, is commonly estimated using either the tracer binding level in a reference region that is devoid of the target [2,3], or a blocking study which involves a baseline PET scan and a second scan with a blocking drug administered just before the tracer [4,5].
In a valid reference region, the tracer is either free or only "non-specifically" bound, and its volume of distribution (V T ) in such a region (V T-RR ) is typically assumed to represent the brain-wide V ND .For many targets this approach is not appropriate because there is no valid reference region, as the target is present throughout the brain [6][7][8][9][10][11][12][13][14][15][16], and thus the signal in any region includes some specific binding.Using an invalid reference region over-estimates V ND , causing underestimation of binding potentials [16], and can confound interpretation of results [17][18][19][20].Automatic extraction of a reference region signal using cluster analysis [7,10,12] of the brain PET data is often not successful, or greatly depends on the data used to train the clustering algorithm [21].
Alternatively, a blocking study with tracer injections before and after a saturating dose of an antagonist with high affinity for the same target of interest allows estimation of brain-wide V ND using a Lassen plot [4,5].However, performing a blocking study in each subject is cumbersome, costly, doubles the radiation exposure, can involve side effects related to the blocking agent, and is therefore generally avoided in clinical research.
A parametric pseudo-reference tissue model was proposed [22] for tracers that have no ideal reference region, which provides estimates only for the binding potential BP ND [2] and not for V ND and thus not for binding potentials BP P and BP F [2], and assumes that BP ND in the pseudo reference region can be estimated from additional competition data.A genomic plot was also recently proposed, which provides V ND estimates only at the population level and requires that the brain maps of messenger RNA transcripts of the specific target of interest be available from the Allen Brain Atlas [23].
Based on compartment models (CMs) [24], we proposed to perform at the individual subject level simultaneous estimation of a common V ND across regions [25] when no valid reference region is available.However, for some tracers such as [ 11 C]DASB (target: serotonin transporter), the simultaneous estimation of V ND across regions often fails to give a unique solution.Separately, we also showed [26] that nonparametric deconvolution is an alternative quantification approach for PET data, which computes binding potentials comparable to estimates by CMs, and for some tracers, shows superior test-retest performance than quantification by CMs [26].
We now propose a new hybrid deconvolution approach (HYDECA) that combines deconvolution and simultaneous search across regions to calculate a brain-wide V ND when arterial blood data are available but a valid reference region is not.HYDECA is validated for [ 11 C] DASB and [ 11 C]CUMI-101 (target: serotonin 1A receptor) using simulations and blocking studies [11,27], and evaluated in test-retest datasets [10,28].

Human subjects and animal studies
Data from published blocking studies in baboons [27] and humans [11], and test-retest datasets in humans [10,28] were used.Human studies were performed in accordance with the 1964 Declaration of Helsinki and its later amendments and approved by The Institutional Review Boards of Columbia University Medical Center (CUMC) and New York State Psychiatric Institute (NYSPI).Animal studies were performed with the approval of the CUMC and NYSPI Institutional Animal Care and Use Committees, according to all applicable regulations governing the use of animals in research.

Nonparametric quantification
According to the extended indicator dilution theory [26,29], the tracer signal in tissue in a brain region i, C Ti (t), after correction for the presence of tracer in vasculature, is a scaled convolution between the metabolite-corrected input function in the arterial plasma, C P (t), and the so-called tissue residue function, R i (t): While K i [mLÁcm -3 Ámin -1 ] is a proportionality constant, R i (t) is defined in the theory of the indicator-dilution method as the fraction of indicator that remains in the tissue after an idealized bolus input concentration at time zero.Initially, the residue must be unity (R i (0) = 1) and from there it decreases (or at least does not increase) with time (refer to [29] for details).Among many nonparametric approaches that can be used to estimate the impulse response function (IRF) in each region i, IRF i (t) = K i R i (t), from known C P (t) and C Ti (t), we proposed using singular value decomposition (SVD) with data-driven selection of the threshold that we described elsewhere [26].

Hybrid deconvolution approach
In the context of PET reversible radiotracers, R i (t) can be interpreted as the fraction of tracer molecules remaining in the tissue over time, and these molecules can be specifically bound to the target, free in water or bound to other molecules.HYDECA decomposes each region R i (t) into the sum of a parametric non-displaceable component, which is approximated as a monoexponential function depending on V ND , assumed common across regions (see details below and comments on the validity of this approximation in the Discussion), and a nonparametric specific component.For any choice of V ND defining the non-displaceable component, the nonparametric specific component can be estimated by subtraction.
Performing such a decomposition for observed PET data can be challenging, but the goal of HYDECA is to objectively ascertain a "reasonable" V ND value by examining data across regions.To illustrate this idea, Fig 1 shows R i (t) curves in two representative regions, calculated using a two-tissue CM (2TCM) [24] (see Eq 8) and based on kinetic rates derived from data with [ 11 C] DASB [28].Non-displaceable component curves based on two "unreasonable" choices of V ND (a value that is 1/4 and 4 times the magnitude of the true V ND , respectively) are compared to the non-displaceable component calculated with the true V ND (the "most reasonable" choice).The non-displaceable component and the corresponding specific component differentially contribute to R i (t) and two effects can be observed.The first effect is that, at time zero, the difference between the slope of R i (t) and that of the non-displaceable component is small, if the V ND value is close to the true V ND .The second effect is that, when the V ND value used for the nondisplaceable component is larger than the true V ND , the corresponding specific component results in negative values, violating its positivity constraints.Different regions show different contributions of non-displaceable and specific component to R i (t).HYDECA is based on finding a V ND value that, across regions, provides the best compromise between these two effects.
To do so, HYDECA requires as input C P (t) and C Ti (t) curves from a pre-determined set of N brain regions, once corrected for the presence of vasculature.In our implementation we assumed a brain-wide blood volume of 5%.HYDECA estimates V ND as follows: IRF i (t) is estimated in each region i from C Ti (t) and C P (t) using SVD as described [26]; K i is obtained as the value of IRF i (t) at time zero (R i (0) = 1 by definition for an idealized bolus input; see "Implementation" section below for comments); R i (t) is then obtained dividing IRF i (t) by the K i estimate; R i (t) is expressed in each region i as the sum of a parametric non-displaceable component (corresponding to an ideal one-tissue CM with distribution volume of V ND ), R ND (t), and a nonparametric specific component, S i (t): Assuming a mono-exponential for R ND (t) represents an approximation (if a 2TCM is needed to describe the data in a given region, R ND (t) would be described by two exponentials [24]), whose validity varies across regions (see comments on the validity of this approximation in the Discussion); HYDECA examines data across regions to find a suitable common V ND .The property expressed in Eq (2) can be derived from CMs as follows: with binding potential BP P as in [2]; HYDECA expresses parametrically only R ND (t).
For fixed values of the tuning parameters β and γ, the following cost function is minimized over V ND using all N regions: Minimization of the first term in Eq (6), which represents the residual sum of squares between R i (t) and R ND (t), calculated up to time γ after tracer injection, identifies V ND values that provide R ND (t) curves with a slope close to the slope of R i (t) at time zero.Given the difficulty in accurately estimating slopes from noisy data, the difference in the slope is approximated as the residual sum of squares between the two curves.The tuning parameter γ controls the number of data points considered for this calculation.Minimization of the second term, which represents the negative area of the curve of the corresponding S i (t), in the case that a portion of S i (t) assumes negative values, penalizes V ND values that lead to unphysiological S i (t) values (Fig 1).If S i (t) is everywhere positive then there is no contribution of the second term.If S i (t) has negative values, then the time t Ã is data-derived as the time point after which S i (t) has consistently positive values.The tuning parameter β weights the contribution of the second term relatively to the first term.We propose and compare two strategies for setting optimal values for the tuning parameters β and γ for a given tracer.

Tuning with simulations
One strategy involves simulating data that imitate characteristics of real data for the tracer at hand, letting β and γ vary over a grid of possible values, and identifying optimal β and γ as  those values that allow HYDECA to generate an estimate of V ND that, on average across all simulated instances, is closest to the true simulated V ND (V ND_TRUE ).We considered a metabolite-corrected input function C P (t) and kinetic rate values in the same brain regions we considered in previous publications [25,26] based on available data [10,28]: cerebellar gray matter (CGM), temporal lobe (TEM), hippocampus (HIP), dorsal caudate (DCA), amygdala (AMY), and ventral striatum (VST), for [ 11 C]DASB; CGM, HIP, TEM and occipital lobe (OCC), and cingulate (CIN) for [ 11 C]CUMI-101.Noise-free C Ti (t) curves were generated for each region using a 2TCM [3,24]: where K 1i , k 2i , k 3i , and k 4i are the values for the kinetic rate parameters of region i.Table 1 lists the kinetic rate values used in each of two simulated cases per tracer: 1) common V ND_TRUE is 3, and 50% of the tracer V T-RR is specific binding (cerebellar grey matter V T ~6); 2) common V ND_TRUE is 5, and ~17% of the tracer V T-RR is specific binding (cerebellar grey matter V T ~6).In all cases, we simulated Gaussian noise with zero mean.To ensure realistic noise characteristics, the variance-covariance matrix used to generate simulated noise was estimated from a matrix of residuals, standardized across time points, from the fits for the considered regions using available data [10,28].In all cases, we simulated 1000 C Ti (t) curves for each region.
For each tracer and V ND_TRUE case, we then: 1) considered a grid of β (0.5 to 14; step: 0.5) and γ values (1 to 30 minutes after tracer injection; step: 1); 2) calculated the cost function (Eq 6) corresponding to all combinations of β and γ within the grids, and over a grid of V ND values (0.1 to 7; step: 0.1), in each of the simulated instances; 3) considered the average cost function (across instances) corresponding to each of the combinations of β and γ; 4) estimated V ND as the value that minimizes each of these average cost functions; and 5) calculated the corresponding absolute estimation error as |V ND_TRUE −V ND |.After obtaining the association between each combination of β and γ within the grids and the corresponding bias of the V ND estimate, we selected as optimal β and γ derived via simulations (β opt-S , γ opt-S ) the values providing the smallest bias.

Tuning with blocking studies
Another strategy involves using blocking studies, if available, letting β and γ vary over a grid of possible values, and identifying optimal β and γ as those values in correspondence of which HYDECA provides a V ND that is, on average across all subjects in the dataset, the closest to the V ND estimated using both scans before and after blocking and Lassen plot [5] (V ND_LASSEN ).
We examined 13 healthy controls imaged with [ 11 C]DASB before and after administration of sertraline [11], and 8 pairs of scans performed on healthy baboons with [ 11 C]CUMI-101 before and after either WAY100635 or 8-OH-DPAT [27].
In each pair, we computed V ND_LASSEN using both scans before and after blocking and the same regions considered in simulation.We then: 1) considered the same grids for β and γ as in tuning with simulations; 2) calculated the HYDECA cost function corresponding to all combinations of β and γ within the grids, and over a grid of V ND values (0.1 to 30; step: 0.1), using in each pair only the scan before blocking and the same regions considered in simulation; 3) estimated V ND as the value that minimizes each of these cost functions; 4) calculated the corresponding absolute estimation error as |V ND_LASSEN −V ND |; and 5) calculated the average (across all subjects within a tracer) estimation error obtained for each combination of β and γ.After obtaining the association between each combination of β and γ within the grids and the corresponding bias in the V ND estimate, we selected as optimal β and γ derived via blocking studies (β opt-B , γ opt-B ) the values providing the smallest bias.
In each scan before blocking, we calculated the percent difference (PD VND ) between V ND estimated using HYDECA (with tuning parameters set with either strategy) and the corresponding V ND_LASSEN , as PD VND = 100Á|V ND_LASSEN −V ND |/ V ND_LASSEN .

Implementation
HYDECA, implemented in Matlab R2012b (www.mathworks.com/), is a fast algorithm that runs in ~14 seconds for one subject on an iMac machine, 3.5 GHz Intel Core i7 Processor, once β and γ are determined.The most computationally demanding component is the data- driven selection of SVD threshold [26].The computational time required to optimize the tuning parameters initially for a given tracer depends on the selected strategy.If this is done using simulations, this can take up to a few hours.Using blocking studies, the computation is complete within a few minutes.
Only with an idealized bolus input does R i (t) reach its maximum at time zero, and in such a case, K i could be derived from the value of the reconstructed IRF i (t) = K i R i (t) at time zero.With a realistic bolus infusion of the tracer, R i (t) reaches its maximum at some time t > 0, and implementation that estimates K i as the maximum of the reconstructed IRF i (t) is preferable.Furthermore, in our implementation, all deconvolved R i (t) curves are first shifted to have their maximum value correspond to time zero before calculating the HYDECA cost function in Eq (6).We do not perform any correction for a physiological delay between C Ti (t) in the different regions and C P (t).

Estimation using the non-ideal reference region
To investigate the bias of HYDECA V ND estimates relative to estimates measured using the Lassen plot, and in comparison to the common practice of setting V ND equal to V T-RR even when the reference region is known not to be valid, we utilized only the scans before blocking in the two available datasets to calculate V T in CGM starting from C P (t) and C Ti (t), using both a 2TCM [24] and Likelihood Estimation in Graphical Analysis (LEGA) [30].CGM was chosen as reference region as it has the lowest V T [10,28] and least displacement of all regions examined in our blocking studies [11,27].LEGA provides the best test-retest reproducibility over analysis with CMs and other graphical approaches for estimates with both tracers [10,28].PD VND with respect to V ND_LASSEN was also calculated for V T-RR obtained with both 2TCM (V T-RR,2TCM ) and LEGA (V T-RR,LEGA ).

Application to test-retest data
As V ND is estimated in order to calculate binding potentials, we considered two available testretest datasets with [ 11 C]DASB [28] and [ 11 C]CUMI-101 [10] and investigated the reproducibility of binding potentials derived using HYDECA versus using the purported reference region (CGM).Both test-retest datasets included only healthy controls, who were imaged with the radiotracer in question twice in one day (once in the morning, once in the afternoon) in a test-retest study design.In all scans, we calculated V ND (using HYDECA with optimal β and γ set with either strategy, and considering the same regions used in simulation), V T-RR,2TCM , and V T-RR,LEGA .For each test-retest pair and region, we calculated the percent difference PD VND-TRT as 100 jV TEST À V RETEST j ðV TEST þV RETEST Þ=2 , where V TEST is the V ND or V T-RR estimate in the test scan, and V RETEST the V ND or V T-RR estimate in the retest scan.We compared PD VND-TRT values obtained from the different methods using a two-tailed paired t-test, considering all possible pairwise combinations of methods.
For each test-retest pair and region, we calculated the percent difference for the binding potentials (PD BPP and PD BPND ) as 100 jBP T À BP RT j ðBP T þBP RT Þ=2 , (BP T : test estimate; BP RT : re-test estimate), and computed average and standard deviation (SD) (across subjects within a tracer) of PD BPP and PD BPND values in each region.We compared PD BPP and PD BPND values obtained from the different methods using a two-tailed paired t-test, region by region, considering all possible pairwise combinations of methods.

Simulation studies
Tuning parameters optimization.Optimization of β and γ using simulations and effects of β and γ values on V ND estimates obtained by HYDECA are shown in Fig 2 .As γ (number of points considered for the first term in Eq (6)) increases, β needs to correspondingly increase to weight more the second term in Eq (6)), in order to minimize bias in V ND estimation.
For each tracer and V ND_TRUE case, we selected the optimal β and γ values in correspondence of which HYDECA provides the least biased estimation of V ND (Fig 2 When only 17% of the V T in the non-ideal reference region is specific binding (V ND_TRUE = 5), for both tracers the number of combinations of β and γ in correspondence of which HYDECA provides an average estimation error smaller than using V T-RR is reduced.However, HYDECA with optimized β and γ still generated a robustly more accurate estimate of V ND than V T-RR ([ 11  Cost functions and estimation bias with optimized tuning parameters.HYDECA cost function curves (Eq 6) using (β opt-S , γ opt-S ) as determined using simulations are convex and unimodal (Fig 3).The corresponding distributions of V ND estimates show a bias, calculated as the average of (V ND_TRUE −V ND ) across instances, of -0.008 (V ND_TRUE = 3) and -0.024 (V ND_TRUE = 5) ([ 11    Estimation bias with optimized tuning parameters.Application of HYDECA to individual scans in the blocking studies, with β and γ optimized using either strategies, in comparison to the use of V T-RR is shown in Fig 5 .V ND estimates by HYDECA with either sets of tuning parameters are considerably less biased, relative to V ND estimates from the Lassen plot, than those using V T-RR .Estimation of V ND using 2TCM in the non-ideal reference region is more biased than that obtained by HYDECA and LEGA for both considered tracers.Average (± SD) PD VND values across subjects are: 15.48% (± 9.82) using HYDECA with (β opt-B , γ opt-B ), 15.40% (± 11.65) using HYDECA with (β opt-S , γ opt-S ), 44.16% (± 22.52) using V T-RR,LEGA , and 70.04% (± 24.00) using V T-RR,2TCM ([ 11 C]DASB); 27.81% (± 19.03) using HYDECA with (β opt-B , γ opt-B ), 26.08% (± 17.24) using HYDECA with (β opt-S , γ opt-S ), 70.26% (± 42.82) using V T-RR,LEGA , and 76.10% (± 56.03) using V T-RR,2TCM ([ 11 C]CUMI-101).All V ND and V T-RR estimates for all approaches and both blocking datasets are reported in Table 2.For both tracers, average (across subjects within each tracer) V ND estimates by HYDECA, with β and γ optimized using either strategies, are closer than both LEGA and 2TCM to average values calculated using Lassen plot, which is considered standard in the field for in vivo estimation of V ND , and SD values are overall lower than those for LEGA and 2TCM.

Test-retest studies
Average (across subjects within each tracer) estimates of V ND and V T-RR in the test-retest datasets (Table 3) are consistent with corresponding values found in the blocking datasets (Table 2), although in the case of [ 11 C]CUMI-101 the two datasets are in different species.V ND values we obtain with HYDECA average 27% of total binding in ventral striatum for [ 11 C] DASB, and 22% of the total binding in hippocampus for [ 11 C]CUMI-101, which is generally in line with reports for other PET tracers [31,32].
Test-retest PD VND-TRT values for V ND (Table 4) from the different methods are not statistically significantly different from each other, with the exception of [ 11 C]DASB, in which case PD VND-TRT values obtained by HYDECA are statistically significantly higher (indicating worse reproducibility) than those derived by LEGA (p-values: 0.003 with β and γ set via simulation; 0.002 with β and γ set via blocking study).See Discussion for factors affecting the reproducibility of V ND by HYDECA.
Reproducibility of the binding potentials estimated using HYDECA V ND , with β and γ optimized using either strategies, is compared to that of binding potentials based on V T-RR,LEGA , V T-RR,2TCM , or direct estimation by 2TCM in Fig 6 .PD BPP values obtained using HYDECA with either sets of optimized tuning parameters are close to each other and comparable to values obtained using V T-RR,LEGA .PD BPP values from the different methods are not statistically significantly different from each other, with the exception of 2TCM direct estimation, where PD BPP values are statistically significantly higher (indicating worse reproducibility) than those of all other methods in the case of [ 11 C]DASB in all brain regions except HIP (range of p-combinations of β and γ for which HYDECA provides an average absolute error in the estimation of V ND that is higher than the absolute error committed by using the V T in the CGM as an estimation of V ND (|err RR |), and the ratio between |err RR | and |err HYDb | in correspondence of the optimal β and γ is reported.The yellow and pink circle indicates the optimal combination of the tuning parameters (β opt-S , γ opt-S ) derived using simulation with V ND_TRUE = 3 and V ND_TRUE = 5, respectively, and the ratio between |err HYDs | (in correspondence of β opt-S and γ opt-S ) and |err HYDb | (in correspondence of β opt-B and γ opt-B ) is reported.V ND : non-displaceable distribution volume; V ND_LASSEN : V ND estimated using both scans before and after blocking and Lassen plot; V T : tracer total volume of distribution; CGM: cerebellum grey matter.https://doi.org/10.1371/journal.pone.0176636.g004Fig 5 .Estimation bias with tuning parameters β and γ optimized using blocking studies.Difference between V ND_LASSEN and V ND estimates obtained by HYDECA with (β opt-B , γ opt-B ) set using blocking studies (y-axis; first row), between V ND_LASSEN and V ND estimates obtained by HYDECA with (β opt-S , γ opt-S ) set using simulations (y-axis; second row), between V ND_LASSEN and V ND estimated as the V T in the CGM using 2TCM (y-axis; third row), and between V ND_LASSEN and V ND estimated as the V T in the CGM using LEGA (y-axis; bottom), as a function of V ND_LASSEN (x-axis) in individual scans in the blocking study with [ 11 C]DASB (left) and [ 11 C]CUMI-101 (right).Solid lines indicate the average error; dotted lines indicate average error ± 1.96 standard deviation.The zero line is the dotted black line.V ND : non-displaceable distribution volume; V ND_LASSEN : V ND estimated using both scans before and after blocking and Lassen plot; V ND (HYDECA): V ND estimated using HYDECA; V T-RR,2TCM : distribution volume in the non-ideal reference region calculated using 2TCM; V T-RR,LEGA : distribution volume in the non-ideal reference region calculated using LEGA; V T : tracer total volume of distribution; CGM: cerebellum grey matter; 2TCM: twotissue compartment model; LEGA: Likelihood Estimation in Graphical Analysis.https://doi.org/10.1371/journal.pone.0176636.g005

Discussion
HYDECA is a data-driven approach that estimates V ND for each individual based on his/her PET data from multiple brain regions.HYDECA is intended for tracers and targets for which a valid reference region does not exist.If a valid reference region does in fact exist, then binding potentials based on V T-RR or on reference region approaches are likely to be more accurate than those based on HYDECA.

Tuning parameters
HYDECA implementation requires determination of two tuning parameters, herein denoted β and γ, and we propose two possible strategies to make this choice for a given tracer: using data from blocking studies, or realistic simulations of the tracer in question.It should be noted that using the same tuning parameters across subjects imaged with the same tracer does not result in estimating the same V ND in each subject.
Of the two strategies, the one using blocking studies involves less subjective judgment.When blocking study data are not available, simulations can be used, but simulated V ND values, kinetic rates, and measurement errors should be chosen carefully to obtain realistic representation of the data with the tracer in question.For established tracers, simulations can be set up using kinetic rate values derived from available data or from the literature.For a new tracer, both simulations and validation with blocking studies are recommended.
For [ 11 C]DASB and [ 11 C]CUMI-101, our results suggest that, regardless of the strategy used to optimize the tuning parameters, HYDECA estimates of V ND are considerably less biased than those obtained based on V T-RR .Even with a "sub-optimal" choice of the tuning parameters, HYDECA estimates of V ND are generally less biased than using a non-ideal reference region (Figs 2 and 4).Although the selection strategies can provide different values for β and γ, the resulting bias in the estimation of V ND is similar (|errHYDs|/|errHYDb| ratios in Fig 4).
If we were to optimize β and γ individually for each subject in the blocking datasets, we would observe quite large inter-subject variability in the optimal β and γ: β = 8.85   in the dataset that is used for tuning parameter selection, using individually optimized β and γ instead of values optimized on average across subjects, as we suggest, would lead to an even less biased V ND estimation.However, the question of which β and γ values to use when applying HYDECA to a subject imaged with the same tracer, but for which a blocking scan is not available, would remain.Individually optimized β and γ values are not obtainable in standard practice.Utility HYDECA is a workable algorithm that can be applied to estimate individual V ND in absence of a reference region or individual blocking data, and could therefore be extremely useful in both clinical and research settings.If the target selected for a given PET application lacks a valid reference region, there is no way to accurately estimate V ND (and consequently specific binding to the target), unless one performs a blocking scan for each subject.HYDECA can provide an alternative convenient quantification approach.For tracers for which HYDECA tuning parameters have already been determined, the published optimized tuning parameters can be used.Otherwise, published blocking studies for the tracer in question would constitute the basis to either tune HYDECA directly (if data are accessible) or to set up a simulation.

Reproducibility
HYDECA estimates of V ND (with either strategy to set tuning parameters) lead to binding potentials estimates with test-retest reproducibility that are comparable to estimates based on V T-RR .Note that average PD BPP values are overall lower when based on HYDECA compared with values based on V T-RR , and not merely because V ND estimates by HYDECA are consistently lower than corresponding V T-RR .Detailed related information is provided in the Supplementary Materials (S3 and S4 Figs, S2 Text).We observe on average worse reproducibility of the estimates based on HYDECA when considering BP ND compared to BP P estimates.Because of the nature of the outcome measure and performance metric used here, BP ND values and their corresponding test-retest performance are more sensitive than BP P to values and changes (in between test and retest scan) in the V ND , which appears at the denominator in the indirect definition of BP ND .When using the V T from an invalid reference region to estimate V ND , reproducibility of that measure depends on, among other factors, how much the tissue time activity curve from that region changes between the test and the retest scans.HYDECA, instead, uses tissue time activity curves from multiple regions to determine V ND , and therefore its test-retest performance is affected, among other factors, by how much the tissue time activity curves from all of these regions change between the test and the retest scans.The test-retest percent difference values for HYDECA V ND (Table 4) are on average worse than those for V T-RR calculated using 2TCM and LEGA, especially in the case of [ 11 C]CUMI-101.Reproducibility performance should be considered when deciding which approach to use in longitudinal studies, while the bias of the approach is more important in group comparisons and cross-sectional studies.

Alternative strategies
If blocking scans are available for a certain tracer, they could be used to estimate a populationbased α = V ND_LASSEN /V T-RR ratio, which could then be used for studies with the same tracer to scale each subject V T-RR in the non-ideal reference region to estimate V ND .We applied such approach to the two available blocking datasets.We found the following V ND_LASSEN / V T-RR,LEGA average (± SD) α ratios: 0.710 (± 0.114) for HYDECA with β and γ set via simulation (p-value: 0.050, 0.024, and 0.036, respectively).Also see comments on BP ND reproducibility in "Reproducibility" section above.From a compartment modeling point of view, however, if there is specific binding in the non-ideal reference region, this would correspond to an additional compartment, which would require a subtraction (rather than a multiplicative adjustment) from the total V T in the region, in order to be properly accounted for.A population-based distance d = V T-RR −V ND_LASSEN can be derived if blocking scans are available for a certain tracer as in the case of the scaled V T-RR .In the two available blocking studies, we found the following average V T-RR,LEGA −V ND_LASSEN distance (± SD) d: 3.04 (± 1.56) for [ 11 C]DASB, and 2.14 (± 0.92) for [ 11 C]CUMI-101.We applied such average distance values to the subjects in the available test-retest datasets to A fixed population-based ratio or distance approach, unlike HYDECA, would not take advantage of the information relative to V ND that is implicitly contained in each individual's PET tissue data across brain regions.Such an approach would rely on blocking studies more heavily than HYDECA, for which tuning parameter selection can alternatively be achieved using simulations.
In the Supplementary Materials (S1 and S2 Figs, S1 Text) we report results obtained on alternative nonparametric binding potentials [26] that can be calculated based on HYDECA, including their test-retest reproducibility and the comparison to 2TCM, LEGA, and alternative strategies.

Choice of regions
Regions that are simultaneously considered should be carefully chosen in all approaches that either take advantage of simultaneous estimation across regions [25,[33][34][35][36][37], or jointly estimate common parameters of interest across regions, like occupancy and V ND in the Lassen plot.For simultaneous estimation approaches to perform well, the regions that are considered should in general have kinetic behavior as distinct as possible [36].Including regions with similar kinetic behavior would serve only to increase the dimensionality of the objective function without adding much useful information [36].The variety in kinetic behavior depends greatly on the tracer at hand.In our previous experience with simultaneous-type estimation with [ 11 C]DASB [33] and [ 11 C]CUMI-101 [25], we had carefully selected regions to represent a broad range of kinetic behavior, while avoiding regions that tend to be noisy.We had also previously assessed the properties of nonparametric quantification in these regions using both simulated and clinical data [26].We are therefore using the same regions in this study.

Choice of deconvolution approach
We used here SVD for its speed and ease of implementation, and have characterized its performance in terms of reproducibility and sensitivity to noise in an earlier publication [26].SVD can however be sensitive to potential delay and dispersion of the injected bolus [38,39].More robust approaches to nonparametric deconvolution [39,40] or functional principal components analysis [41] may further improve HYDECA performance.Here we provide a framework for HYDECA and comparison between different implementations of the algorithm is beyond our scope.

Limitations
Vascular correction.The tracer signal in the brain tissue can be modeled as in Eq (1) only after correction for intravascular activity.Here, following a practice common in the field, we assumed a brain-wide fractional blood volume (V B ) of 5%.It is recommended that the V B value be optimized before applying HYDECA (or any other PET quantification approach) if pathological changes in the fractional blood volume are suspected in the population at hand.HYDECA performance, as that of any PET quantification approach, may in fact be affected by an erroneous choice of the V B value used to correct the tissue time activity curves.We ran an additional simulation to investigate the sensitivity of HYDECA estimates of V ND to a potentially erroneous vascular correction of the measured time activity curves (details are reported in the Supplementary Materials, S5-S7 Figs, S3 Text).HYDECA estimates of V ND appear to be robust to erroneous correction of the time activity curves for errors in V B in the range -4% to +5% for [ 11 C]DASB, and -4% to +7% for [ 11 C]CUMI-101.
As blood volume may vary in the brain, using a brain-wide value may not significantly impact outcome measures such as V T and binding potentials, but may impact the upslope of the tissue signal, and thus the R i (t) estimated nonparametrically.If V B varies across regions, a case that is not trivial for any of the quantification approaches used in PET, one potential strategy to account for this within HYDECA could be incorporating the vascular correction component into the impulse response function that is nonparametrically deconvolved in each region.The problem may be treatable from a mathematical point of view, but would require careful comparison of more sophisticated approaches to deconvolution than SVD.Another potential strategy could be exploiting the semiparametric nature of HYDECA and adding V B as a free-parameter to be estimated in each of the regions that are simultaneously considered, but this would require a more complex optimization procedure than the simple grid approach that we proposed for V ND .Correction for intravascular activity represents just as much of a problem for other approaches proposed as alternatives to compartment models [42,43].
Assumption of a mono-exponential R ND (t) curve.To ensure identifiability of the two components of the residue function curve R(t) (non-displaceable and specific), HYDECA needs to assume a certain shape to describe the non-displaceable component, R ND (t).We chose, in part for its simplicity, a mono-exponential function, which would represent the impulse response function in the case of an "ideal" reference region with total distribution volume equal to V ND .Assuming a mono-exponential curve for R ND (t) represents an approximation: if a 2TCM is needed to describe the data in a given region, the R ND (t) curve of the region would be more appropriately described by a two-exponential function (24).We note that a similar assumption is central in the development of the very widely used simplified reference tissue model (SRTM) [44], which assumes that the total (non-displaceable plus specific) impulse response function of the target region (which, as well, would be a two-exponential function) can be reasonably approximated by a mono-exponential curve.The Supplementary Materials (S8 and S9 Figs, S4 Text) report data to evaluate the validity of such approximation for the two tracers considered here.Our evaluation indicates that a mono-exponential approximation for R ND (t) would be problematic only in the situation in which k 3 >> k 4 , which means that more tracer molecules transit in a given amount of time from the non-displaceable binding state into the specific binding state than vice versa.We recommend that the simplifying assumption of a mono-exponential R ND (t) curve be evaluated for tracers for which it is suspected that k 3 >> k 4 .However, we remind the reader: 1) that HYDECA uses data across many regions, for some of which the mono-exponential assumption may hold better than for others, and provides a brain-wide value of V ND that satisfies certain constraints (via the HYDECA cost function) on average across such regions; and 2) that parts of the R ND (t) curve that are potentially erroneously determined in a region due to the simplifying mono-exponential assumption are likely to be captured by the corresponding nonparametric R S (t) curve, for which there is no assumption besides being positive and monotonic.We want also to stress that the assumption of a common, brain-wide V ND implies that the ratio of the transfer constants (V ND = K 1 /k 2 ) is the same everywhere in the brain for non-specific binding.This same assumption is routinely made when using CMs and/or graphical approaches in a reference region to estimate a brain-wide common V ND , when constraining the K 1 and k 2 parameters in a 2TCM to those of a reference region, or when using SRTM.
Applicability to other populations.The two assumptions required to apply HYDECA are that: a) the non-displaceable distribution volume V ND is uniform brain-wide within each subject (which is the same assumption regularly considered in the field when estimating V ND from a reference region, or when using SRTM); and b) the non-displaceable component of the residue function, R ND (t), is reasonably described by a mono-exponential function (a similar assumption is considered for both reference and target region when using SRTM).So unless there is a population or group of subjects where it is suspected that these two assumptions are seriously violated, HYDECA can be applied.The presence of altered kinetics in the tissue time activity curves of such a population would be problematic for any of the other PET quantification approaches that are based on the assumptions above.

Future investigations
Future investigations include developing a method to provide a measure of precision [45] for HYDECA V ND estimates, validating HYDECA across tracers, and assessing whether performing the tuning of β and γ only once for a given tracer will suffice, which should be the case if the noise characteristics and kinetics range of independent data acquired with a tracer for which the HYDECA tuning parameters have been determined will resemble those of the data used in such determination.

Conclusions
We showed, using two PET radiotracers that, in the absence of a valid reference region, HYDECA can provide individual estimates of a brain-wide V ND without requiring a blocking study, and these estimates are less biased, with respect to estimation with Lassen plot, which is the method of reference, as it represents a standard in the field for in vivo estimation of V ND in humans, than those obtained relying on the V T in a non-ideal reference region.

Fig 1 .
Fig 1. Illustration of the idea behind the algorithm in the hybrid deconvolution approach.R i (t) curves (black lines) calculated using the IRF of the 2TCM and values of the kinetic rates derived from a study with [ 11 C]DASB for 2 representative regions.Red solid lines indicate the non-displaceable component calculated with a V ND that is 1/4 the size of (top) and 4 times higher than (bottom) the true underlying V ND ; green dotted lines indicate the nondisplaceable component calculated with the true underlying V ND ; blue lines indicate the corresponding specific component.Inset plots are added to allow closer inspection of the first 10 minutes after tracer injection.The time t* is derived from the data as the time point after which the specific component has consistently positive values.IRF: impulse response function; 2TCM: two-tissue compartment model; V ND : non-displaceable distribution volume.

Fig 2 .Fig 3 .
Fig 2.Optimization of tuning parameters β and γ using simulations.Average absolute error in the estimation of V ND as a function of the values for the tuning parameters β and γ, for all simulated cases (V ND_TRUE = 3, top; V ND_TRUE = 5, bottom) and tracers.Each point in the matrices correspond to a specific combination of β (vertical axis) and γ (horizontal axis) values in the selected grids, and represents the average (across simulated instances) absolute distance between V ND estimate obtained by HYDECA, using the corresponding combination of β and γ, and V ND_TRUE.The white circle indicates the optimal combination of the tuning parameters (β opt-S , γ opt-S ) derived using simulations in each case, and the average absolute error in the estimation of V ND in correspondence of the optimal tuning parameters is reported (|err HYDs |). Green circles indicate the combinations of β and γ for which HYDECA provides an average absolute error in the estimation of V ND that is higher than the absolute error committed by assuming V ND equal to the V T in the nonideal reference region (CGM) (|err RR |), and the ratio between |err RR | and |err HYDs | in correspondence of the optimal β and γ is reported.V ND : nondisplaceable distribution volume; V ND_TRUE : true simulated V ND ; V T : tracer total volume of distribution; CGM: cerebellum grey matter.https://doi.org/10.1371/journal.pone.0176636.g002

Fig 4 .
Fig 4. Optimization of tuning parameters β and γ using blocking studies.Average absolute error in the estimation of V ND as a function of the values for the tuning parameters β and γ for both tracers.Each point in the matrices correspond to a combination of β and γ values in the selected grids, and represents the average (across scans within the same tracer) absolute distance between the V ND estimated by HYDECA, using the corresponding combination of β and γ, and V ND_LASSEN .White circles indicate the optimal combinations of the tuning parameters (β opt-B , γ opt-B ) derived using the blocking studies, and the average absolute error in the estimation of V ND in correspondence of (β opt-B , γ opt-B ) is reported (|err HYDb |). Green circles indicate the ± 4.64, γ = 19.77± 11.80 ([ 11 C]DASB); β = 9.31 ± 5.04, γ = 18.38 ± 12.37 ([ 11 C]CUMI-101).For subjects

Fig 6 .
Fig 6.Reproducibility of binding potentials estimated using HYDECA, LEGA, and 2TCM.Average plus standard deviation (across test-retest pairs of scans within each tracer) test-retest percent difference PD BPP values calculated in each of the considered region for [ 11 C]DASB (left) and [ 11 C]CUMI-101 (right), using BP P based on V ND from HYDECA, BP P based on V T-RR,LEGA , BP P calculated directly from the 2TCM kinetic rates, and BP P based on V T-RR,2TCM (top).Corresponding values for the test-retest percent difference PD BPND (bottom).Vertical axes are reported in logarithmic scale to allow for easier visualization of the direct 2TCM results based on 2TCM kinetic rates.Statistically significant comparisons (p-value 0.05) are indicated.AMY: amygdala; CIN: cingulate; DCA: dorsal caudate; HIP: hippocampus; OCC: occipital lobe; TEM: temporal lobe; VST: ventral striatum; PD BPP : percent difference for BP P ; PD BPND : percent difference for BP ND ; V ND : non-displaceable distribution volume; V T-RR,LEGA : tracer total distribution volume in the nonideal reference region estimated using LEGA; LEGA: Likelihood Estimation in Graphical Analysis; 2TCM: two-tissue compartment model; V T-RR,2TCM : racer total distribution volume in the non-ideal reference region estimated using 2TCM.https://doi.org/10.1371/journal.pone.0176636.g006 [ 11 C]DASB, and 0.619 (± 0.145) for [ 11 C]CUMI-101.We applied such ratios to the subjects in the available test-retest datasets to calculate a scaled V T-RR,LEGA , and then calculated the corresponding BP P-α = V T (LEGA)-αV T-RR,LEGA and BP ND-α = BP P-α /αV T-RR,LEGA values and their test-retest percent difference (Fig 7).Test-retest percent differences values obtained using the different methods reported in Fig 7 are not statistically significantly different from each other in the case of BP P , nor in the case of BP ND , with the exception of [ 11 C]DASB BP ND in DCA, and [ 11 C]CUMI-101 BP ND in OCC and CIN, for which percent differences values obtained using a population-based α ratio are statistically significantly lower (indicating better reproducibility) than those based on

Fig 7 .
Fig 7. Reproducibility of binding potentials estimated using HYDECA and alternative strategies based on blocking studies.Average plus standard deviation (across test-retest pairs of scans within each tracer) test-retest percent difference PD BPP values calculated in each of the considered region for [ 11 C]DASB (left) and [ 11 C]CUMI-101 (right), using BP P based on V ND from HYDECA, BP P based on scaled V T-RR,LEGA , and BP P based on average distance d (top).Corresponding values for test-retest percent difference PD BPND (bottom).Statistically significant comparisons (p-value 0.05) are indicated.AMY: amygdala; CIN: cingulate; DCA: dorsal caudate; HIP: hippocampus; OCC: occipital lobe; TEM: temporal lobe; VST: ventral striatum; PD BPP : percent difference for BP P ; PD BPND : percent difference for BP ND ; V ND : non-displaceable distribution volume; V T-RR,LEGA : tracer total distribution volume in the non-ideal reference region estimated using LEGA; LEGA: Likelihood Estimation in Graphical Analysis.https://doi.org/10.1371/journal.pone.0176636.g007

S5Fig.
Sensitivity of HYDECA estimates of the non-displaceable distribution volume (V ND ) to erroneous vascular correction.Percent difference (PD errVC ) between the nondisplaceable distribution volume (V ND ) value estimated at each instance of erroneously corrected time activity curves and the V ND value estimated in correspondence of the accurately corrected set of time activity curves (y-axis), as a function of the difference between the true fractional blood volume (V B ) value and the value adopted for correction (x-axis); dots and error bars indicate average and standard deviation across subjects, respectively, within each tracer.The dotted horizontal lines indicate the +10%, 0%, and -10% mark, respectively.(PDF) S6 Fig. Residue function curves R(t) and vascular correction: [ 11 C]CUMI-101.Residue function curves R(t) in correspondence of different errors and no error in the fractional blood volume value (V B ), and the corresponding HYDECA cost functions, in a representative subject for [ 11 C]CUMI-101.CIN: cingulate; HIP: hippocampus; OCC: occipital lobe; TEM: temporal lobe; CGM: cerebellum grey matter.(PDF) S7 Fig. Residue function curves R(t) and vascular correction: [ 11 C]DASB.Residue function curves R(t) in correspondence of different errors and no error in the fractional blood volume value (V B ), and the corresponding HYDECA cost functions, in a representative subject for [ 11 C]DASB.AMY: amygdala; DCA: dorsal caudate; HIP: hippocampus; TEM: temporal lobe; VST: ventral striatum; CGM: cerebellum grey matter.(PDF) S8 Fig. Validity of mono-exponential assumption for the residue function non-displaceable component R ND (t): [ 11 C]DASB.Average (across time points) square distance between the residue function non-displaceable component, R ND (t) (see S4 Text), with k 3 and k 4 >0, and R ND (t) with k 3 = k 4 = 0 as k 3 and k 4 vary, in 4 cases of (K 1 , k 2 ) for [ 11 C]DASB.VST: ventral striatum; CGM: cerebellum grey matter.K 1 , k 2 k 3 and k 4 : kinetic rate parameters of a two-tissue compartment model.(PDF) S9 Fig. Validity of mono-exponential assumption for the residue function non-displaceable component R ND (t): [ 11 C]CUMI-101.Average (across time points) square distance between the residue function non-displaceable component, R ND (t) (see S4 Text), with k 3 and k 4 >0, and R ND (t) with k 3 = k 4 = 0 as k 3 and k 4 vary, in 4 cases of (K 1 , k 2 ) for [ 11 C]CUMI-101.HIP: hippocampus; CGM: cerebellum grey matter.K 1 , k 2 k 3 and k 4 : kinetic rate parameters of a two-tissue compartment model.(PDF) S1 Text.Alternative nonparametric binding potentials and their test-retest reproducibility.Supporting information and equations accompanying S1 and S2 Figs.(PDF) S2 Text.Comparison of binding potentials and test-retest percent difference values.Supporting information and equations accompanying S3 and S4 Figs.(PDF) S3 Text.Sensitivity to vascular correction.Supporting information and equations accompanying S5-S7 Figs.(PDF) S4 Text.Assumption of a mono-exponential non-displaceable residue function.Supporting information and equations accompanying S8 and S9 Figs.(PDF)

Table 1 . List of kinetic rate values used in the simulations. region K 1 [mLÁcm -3 Ámin -1 ] k 2 [min -1 ] k 3 [min -1 ] k 4 [min -1 ] V T [mLÁcm -3 ] [ 11 C]DASB V ND = 3
https://doi.org/10.1371/journal.pone.0176636.t001 Similarly, PD BPND values from the different methods are not statistically significantly different from each other, with the exception of 2TCM direct estimation in the case of[ 11 C]DASB in all regions (PD BPND values statistically significantly higher than those of all other methods; range of p-values: 1.10E-5 to 0.032), and 2TCM direct estimation in the case of[ 11C]CUMI-101 in TEM and CIN, for which PD BPND values are statistically significantly higher than those derived by LEGA (p-values: 0.035 in TEM; 0.021 in CIN) and by 2TCM indirect estimation (pvalue: 0.032 in CIN).Overall, the test-retest reproducibility of binding potentials obtained using all methods reported in Fig 6 are comparable, with the exception of 2TCM direct estimation in the case of [ 11 C]DASB.