Estimation of in vivo nonspecific binding in positron emission tomography studies without requiring a reference region

NeuroImage 108 (2015) 234–242

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Estimation of in vivo nonspecific binding in positron emission tomography studies without requiring a reference region R. Todd Ogden a,b,⁎, Francesca Zanderigo b, Ramin V. Parsey c a b c

Department of Biostatistics, Columbia University, 10032, New York, NY, USA Department of Molecular Imaging and Neuropathology, New York State Psychiatric Institute, 10032, New York, NY, USA Departments of Psychiatry and Radiology, Stony Brook University, Stony Brook, 11794, NY, USA

a r t i c l e

i n f o

Article history: Accepted 14 December 2014 Available online 24 December 2014 Keywords: Kinetic modeling Reference region Simultaneous estimation

a b s t r a c t Estimation of outcome measures in in vivo neuroreceptor mapping with positron emission tomography (PET) commonly depends on an assumption of uniform nondisplaceable binding throughout the brain. In many cases, this can be estimated based on data from a “reference region,” an area of the brain devoid of the receptor of interest. However, often such a region does not exist, as there are some receptors everywhere throughout the brain. Erroneously designating a region as a “reference” can lead to biased estimation, and furthermore, if the level of specific binding in the purported reference region differs between comparison groups, the validity of resulting conclusions may be called into question. We present a method for estimation of all common PET outcome measures that can provide good estimates even when no reference region exists. Our aim is to use information from several regions simultaneously to estimate the information common to all regions. By not requiring specification (or validation) of a reference region, such an approach can provide an automated, objective approach for kinetic modeling of PET data. We illustrate the performance of these methods on simulated data, human [11C]WAY-100635 data, and [11C]CUMI-101 blocking data in baboons. We show close agreement between estimates obtained by using the proposed method (which does not require a reference region) and estimates based on either a reference region or a blocking study. © 2014 Elsevier Inc. All rights reserved.

Introduction In order to more accurately assess diagnosis and risk of various psychiatric and neurological diseases and ultimately to develop more effective therapeutic interventions, it is necessary to improve our understanding of the pathophysiology of these conditions. An invaluable tool for this is positron emission tomography (PET), a nuclear imaging technology that uses radioactively labeled molecules (radioligands) to quantify, for example, the distribution of proteins throughout the brain. In neuroreceptor mapping applications with PET, the density of a receptor of interest throughout the brain can be estimated by modeling the kinetic behavior of a radioligand that binds with the receptor of interest, with density measures being derived from estimated model parameters. One very common model used for this purpose is the two-tissue compartment model (Innis et al., 2007), diagrammed in Fig. 1, in which molecules may move among three compartments. The first compartment is termed the “plasma compartment” in which a

⁎ Corresponding author at: Department of Biostatistics, Columbia University, 6th floor 722 W. 168th St., New York, NY 10032, USA. Tel.: +1 212 342 1247; fax: +1 212 305 9408. E-mail address: [email protected] (R. Todd Ogden).

http://dx.doi.org/10.1016/j.neuroimage.2014.12.038 1053-8119/© 2014 Elsevier Inc. All rights reserved.

molecule moves throughout the body in the bloodstream. To enter the other compartments, a molecule must cross the blood brain barrier (BBB), and when it does it may bind to the target receptor (the “specifically bound” compartment). Nondisplaceable (ND) uptake includes both radioligand nonspecifically bound (bound to molecules other than the target of interest) and free radioligand in the tissue (Innis et al., 2007). During PET imaging it is possible only to measure the total concentration of the radioligand in the brain, representing the total of both the specifically bound and the free/nonspecifically bound compartments. Using blood samples drawn from an artery during scanning, it is possible to measure the concentration in the plasma compartment over time, the fraction of the unmetabolized compound in the blood, and the “free fraction,” the fraction not bound to plasma proteins. The parameters K1, k2, k3, and k4 indicate the rates at which molecules move between compartments. If a region has no receptors of interest, then there is no “specifically bound” compartment in Fig. 1 and the appropriate model reduces to a one-tissue compartment model that depends on only two rate parameters, typically denoted K1′ and k2′ if the corresponding region is to be used as a reference region. Based on these kinetic models, which are based on an assumption of first-order kinetics throughout, and given knowledge of the rate parameters, various macroparameters can be computed, including the total

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

235

Fig. 1. The two-tissue kinetic model.

volume of distribution, which for the two compartment model is defined as VT ¼


 K1 k 1 þ 3 ¼ V ND þ V S : k2 k4

ð1Þ

This total consists of two parts: VS is the “displaceable” or specific volume; and VND is the “nondisplaceable” volume, which corresponds to the free radioligand and the nonspecific binding. In the two-tissue model these components may be computed as VND = K1/k2, and thus the specific volume of distribution is the difference in these two quantities, (K1k 3)/(k 2k 4). In the one-tissue model, there is no displaceable volume, and so the total volume of distribution consists only of the nondisplaceable volume. Measures related to the density of neuroreceptors may be computed from these volumes, including BPP:

fact constant throughout the brain, it is possible to improve estimation of binding measures by identifying a reference region, an area of the brain that is known to be devoid of the receptor of interest. If such a region may be found, as described above, a one-tissue compartmental model should be appropriate for that region, and for such a model the total distribution volume is simply the ratio of K1′ and k2′. Again, assuming that the non-displaceable volume is constant throughout the brain, the parameter estimates from the well-identified model for the reference region may be used to constrain the corresponding parameters in the other regions/voxels. In particular, we could restrict the ratio K1/k2 to be constant for all regions, and in fact, we could constrain this ratio ^ ND , the estimated ratio of these estimated quantities to be equal to V from the reference region. Model fitting for any region/voxel thus requires minimizing Xn i¼1

BP P ¼

K 1 k3 ; k2 k4

ð2Þ

as well as BPF = BPP/fP, where fP represents the plasma free fraction. Defining the parameter vector θ = (K1, k2, k3, k4) the two-tissue kinetic model specifies that the concentration of the radioligand at time t is given by f ðt; θÞ ¼

2 X

  L j ðθÞ exp −R j ðθÞt ⊗C P ðt Þ

j¼1

where CP(t) is the “arterial input function”, i.e., the concentration of the compound in plasma at time t, that has been corrected for the fraction of radioligand that does not penetrate the BBB because it is bound to the blood proteins. The parameters are L1, L2, R1 and R2, which are each functions of the rate parameters K1, k2, k3 and k4. Given plasma data and concentration data from the PET scanner representing a timeactivity curve (TAC) (Y1, t1), …, (Yn, tn) where Y i ¼ f ðt; θÞ þ ε i ;

i ¼ 1; …; n;

ð3Þ

with εi distributed independently having mean zero and Var(εi) = σ 2/wi, the kinetic rate parameters may be estimated using standard nonlinear regression methodology, i.e., calculating ^θ ¼ argmin θ

n X

2

wi ðY i −f ðt i ; θÞÞ

i¼1

Here we suppose the variance factors (w1, …, wn), which take into account frame duration, radioligand decay rate, region size, etc., are known. The individual rate parameters are typically not well identifiable, and in many situations, even macroparameters like BPP in Eq. (2) are not identifiable, making resulting estimation of binding outcome measures unstable. A popular and very general approach to this problem requires the assumption that the level of non-displaceable binding is constant throughout the brain. This assumption can be examined for a given radiotracer using blocking studies. If non-displaceable binding is in

2

wi ðY i − f ðt i ; θÞÞ

ð4Þ

^ . over θ under the constraint that the first element of θ is equal to k2 V ND This reduces from four to three the number of free parameters that must be estimated for each region of interest (ROI)/voxel, and can greatly increase the stability of the estimation of outcome measures (Parsey et al., 2000). (Here we emphasize that our use of the term “reference region” in this development is distinct from the popular “reference tissue models” that seek to estimate BPND when arterial data are not available.) Alternatively, if a reference region is present, allowing direct estimation of VND, “indirect” estimation of binding measures can be accomplished via subtraction, i.e., ^ −V ^ ; c ¼V BP P T ND

ð5Þ

^ is the estimate according to Eq. (1) for the region/voxel comwhere V T ^ is the puted without placing any constraints on the parameters and V ND

estimated total volume of distribution of the reference region as before. We note that these approaches can be combined; i.e., the estimation of VT in Eq. (1) could be done using the earlier constraint. However, the issue of identifying and validating a reference region is problematic. Traditionally, an anatomical reference region is suggested by previous in vitro studies. Given this knowledge, a reference region for each subject may be identified on a higher-resolution magnetic resonance imaging (MRI) scan of each subject and then co-registered to the corresponding PET images. More recently, some attempts, generally either supervised or unsupervised clustering, have been made to identify groups of individual voxels that demonstrate kinetic behavior that resembles what would be expected in a reference region (DeLorenzo et al., 2009; Turkheimer et al., 2007; Yaqub et al., 2012). Setting aside the difficulties of finding and validating a reference region, in many situations, there is some level of displaceable binding throughout the entire brain, and thus no reference region exists. In such a case, if we erroneously take as a “reference region” one that in fact contains some specific binding, resulting estimates of the binding outcome measures will be biased, and the extent of the bias will depend on the level of specific binding in the purported reference region. This would be particularly troublesome when attempting to establish differences in average binding between groups if the level of non-

236

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

displaceable binding differs between groups. And without a reference region, it would be difficult to determine whether this is the case. We propose a straightforward procedure that will allow estimation of all outcome measures without requiring a reference region. The key idea is that information about the level of non-displaceable binding is contained in all TACs, and by considering data from multiple regions at once it is possible to estimate VND without the necessity of designating a reference region. This will allow for a completely automated, objective method for estimating binding parameters. Methods If non-displaceable binding is constant across the brain, then all regions/voxels share a common VND. Rather than fitting each region separately, fitting multiple regions simultaneously allows us to “borrow strength” across all regions, and to combine the pieces of information about VND contained in each region TAC into a stable estimate of this quantity. Simultaneous estimation across regions in kinetic modeling has been applied in other contexts as well, generally to allow more complicated model structures to be fitted identifiably by constraining some parameters to be equal across regions (Zamuner et al., 2012; Huesman and Coxson, 1997; Millet et al., 2006; Raylman et al., 1994). To our knowledge, however, such an approach has never been used to estimate VND. In particular, we propose to constrain K1/k2 = VND to be the same for all regions and simultaneously estimate the three other kinetic parameters specific to each ROI as well as the common parameter VND, i.e., letting θr be the vector of rate parameters for region r and Yir to be the observed concentration in region r and frame i, to estimate VND, θ1,…, θR as the minimizers of XR Xn r¼1

i¼1

2

wir ðY ir − f ðt i ; θi ÞÞ

ð6Þ

under the constraint that the first element of θr is equal to VND times the second element of θr, where the wir values are a set of known weights. (The weights should take into account both the frame duration and the decay of the radioligand as well as the relative noise levels of the various regions, which naturally will depend on the size of the regions. In both our simulation study and our real data applications we take wir to be voxel size divided by squared frame duration (Parsey et al., 2000).) Here, the kinetic parameters for region r are θr = (k2;rVND, k2;r, k3;r, k4;r). This gives a total of 3R + 1 free parameters when fitting R regions: one common VND and three free rate parameters for each region. While this general approach offers a solution for estimating binding outcome measures when no reference region is available, the primary cost is that of greater computational complexity. When a reference region is available and used to constrain the optimization, there is one three-dimensional optimization required for each region. Simultaneous estimation would require only a single optimization that will provide parameter estimates for all R regions, but it would be a much higher dimensional problem than the region-specific optimizations. A variety of general-purpose algorithms are in widespread use for fitting high dimensional models such as Eq. (6). In the situation considered here, there is some structure that can be exploited in order to improve both the likelihood of converging to the global optimum and the speed of the convergence. We note that if the value of VND were given, estimation of the rate parameters for any given region could be accomplished with a standard nonlinear least-squares method in three-dimensional space, as in Eq. (4). In Eq. (6), there is only one parameter (VND) that is common to all regions. Thus optimization of Eq. (6) may be done in a nested fashion, with a one-dimensional optimization at the top level and R 3-dimensional optimizations subordinate to it. In particular, we propose to regard Eq. (6) as a one-dimensional (over VND) problem, noting that, given the current value for VND, calculation of the objective

function involves R separate 3-dimensional optimizations. The nested optimization could be represented as ^ ¼ argmin V ND V ND

nX h R r¼1

minθ

Xn

2

o

w ðY ir −f ðt i ; θr ÞÞ  i¼1 ir

ð7Þ

again with the constraint on the inner minimization that the first element of θr is equal to VND times the second element of θr, noting that the expression in the square brackets is a function of VND. The primary advantage of this approach relative to other potential approaches that we describe in the Discussion is computational stability. It requires only low-dimensional optimizations and may be expected to find the global optimum with high probability. This algorithm is relatively easy to implement, and it is straightforward to parallelize, thereby lending itself well to distributed computing environments. In addition, some understanding of the uncertainty in estimation of VND may be gained by examining a profile plot, which would show possible values of VND on the x-axis and the corresponding optimized criterion (calculated by evaluating the quantity in curly braces in Eq. (7)) on the y-axis. If the minimum of the plot is relatively sharp, that would indicate high identifiability of VND, while if the profile plot is fairly flat, there is more uncertainty in estimating VND (and thus also in the resulting binding outcome measures). Two examples of profile plots are given in Fig. 2. However, since these functions are not always smooth and convex, we prefer instead a one-dimensional grid-search algorithm for optimizing rather than standard optimization algorithms. This approach, while generally less computationally efficient than the nested optimization approach will provide much more reliable estimates in general since it does not rely on the convexity of the profile plot. Thus, provided of course that all parameters are identifiable, implementing a grid-search algorithm as described here can be expected to give reliable estimates in a wide variety of situations. In our applications we consider: for the [11C]WAY data, a grid of 1000 equally spaced VND values ranging from 0.01 to 1; for the [11C]CUMI data (both simulated and on baboons), a grid of 4000 equally spaced VND values ranging from 0.01 to 10; these endpoints having been chosen so as to fully enclose estimates of nondisplaceable binding that arise from analysis of existing data. Results Simulation study To investigate the performance of the two methods, we prepared a simulation study designed to generate data that imitate characteristics of real data. In particular, we used parameters estimated in a study with 11C-CUMI-101, a serotonin 1A receptor partial agonist (Milak et al., 2010, 2011). We selected five regions, the cerebellar grey matter (cgm), the hippocampus (hip), the temporal lobe (tem), the occipital lobe (occ), and the cingulate (cin), chosen to give a range of VT values. Based on existing data on control subjects, we set the VT values for the five regions to be 6 (cgm), 7.8 (occ), 8.4 (cin), 9 (tem), and 16 (hip), respectively. We set the “true” value of VND to be 3. To ensure realistic noise characteristics, we estimated the variance-covariance matrix of the errors at each time point and simulated Gaussian noise having the same covariance structure. We simulated 100 datasets with time-activity curves for each of these five regions and performed estimation of VND as well as the VT values for each region. Results of the simulation are displayed in Fig. 3, and summary statistics are provided in Table 1. The left-hand panel displays smoothed histograms of estimates of VND (black curves) and VT (colored curves) for each of the regions considered. For each smoothed histogram, a vertical line marks the “true” value of the parameter. Estimates of VND and of VT are nearly unbiased. The right-hand panel shows smoothed histograms for estimates of

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

237

Fig. 2. Profile plots for estimating VND using data from two different studies. The vertical lines indicate the respective minima of the two functions.

BPP. With this outcome measure, the histograms of estimates are centered near the corresponding “truth” in each situation. Application to [11C]WAY-100635 human data For demonstrating the performance of the proposed methodology we consider human PET imaging data using the radioligand [11C]WAY-100635 (N-[2-[4-(2-methoxyphenyl)-1-piperazinyl]ethyl]N-(2-pyridyl)cyclohexanecarboxamide), a serotonin 1A receptor (5HT1A) antagonist (Forster et al., 1995). This radioligand is modeled well using a two-tissue model constraining the K1/k2 ratio for each region to be equal to that estimated from the cerebellar white matter (Parsey et al., 2000). For this radioligand, the reference region is well validated (Parsey et al., 2005, 2010; Hirvonen et al., 2007), and this allows for ready comparison between estimates of VND by the proposed method and estimates based on modeling the cerebellar white matter. Weights are calculated based on frame durations (Parsey et al., 2000) and ROI size.

We applied the proposed methodology to data from 44 control subjects imaged with [11C]WAY-100635 with full arterial sampling and imaging protocol as described in Parsey et al. (2010). For application to the proposed method, we selected five regions of interest: the hippocampus, the temporal lobe, the occipital lobe, and anterior cingulate and the thalamus, in order to give a reasonable range of kinetic behavior. Results on estimating VND are shown in Fig. 4. Our estimation method shows strong agreement with the VND estimates obtained using a one-tissue model in the cerebellar white matter. Achieving close agreement on VND estimates is an important goal, but a more relevant consideration is whether the customary kinetic outcome measures can be estimated well when using the proposed methods in the absence of a reference region. Figs. 5 and 6 display the results of estimating VT and BPP, respectively, for several regions obtained by our proposed method for estimating VND, allowing comparison to estimates obtained by exploiting the existence of a reference region. The agreement is generally quite good.

Fig. 3. The left panel displays a smoothed histogram of estimated VND values (black curve) along with estimated VT values for each region considered (color curves). For each smoothed histogram, a vertical line indicating the “truth” is also shown in the same color. The right panels displays smoothed histograms of estimated BPP for each of the five regions, again with a vertical line indicating the “true” value of BPP for each region.

238

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

Table 1 Summary statistics for the simulated data. Volume of distribution

cgm hip tem occ cin Vnd

BPp

Truth

Mean

SD

rootMSE

Skewness

Truth

Mean

SD

rootMSE

Skewness

6.0 16.0 9.0 7.8 8.4 3.0

6.12 16.54 9.17 7.92 8.58 3.04

1.04 1.66 1.04 0.81 1.01 0.66

1.05 1.74 1.05 0.82 1.02 0.66

5.28 1.07 0.10 0.25 0.49 1.09

3.0 13.0 6.0 4.8 5.4

3.08 13.51 6.13 4.88 5.54

1.15 1.30 1.18 0.98 1.13

1.15 1.39 1.18 0.97 1.14

1.94 0.55 −2.56 −2.46 −2.13

Application to [11C]CUMI-101 baboon blocking data We also applied the methodology developed here to data from a [11C]CUMI-101 blocking study on baboons (Milak et al., 2008). As there is displaceable binding in the baboon cerebellum, VND was estimated using our proposed method and compared with that estimated based on the Lassen plot (Lassen et al, 1995). We applied the proposed methodology to data from four scans on two different baboons imaged with full arterial blood sampling and imaging protocol as described in (Milak et al, 2008). In our application, we selected 20 regions of interest. Results on estimating VND are shown in Fig. 7. On the x-axis are estimates made according to the Lassen plot (Lassen et al., 1995) and on the y-axis are the estimates made using our proposed method. Though the number of scans is small, there is generally good agreement in these estimates. Figs. 8 and 9 show the comparison of estimates of VT and BPP, respectively, using the methods described here with those obtained by LEGA (Ogden, 2003; Parsey et al., 2003), a method that been demonstrated to give best results for this radioligand in test-retest studies (Milak et al., 2008, 2010). Again, there is fairly strong agreement throughout. Discussion In routine application, this algorithm (nonlinear least squares for R regions nested within a grid search for VND is necessary only once per scan, and the primary purpose of this is to estimating VND for the scan. Subsequently, VT can be estimated for all regions or voxels using any

b

0.2 0.1 −0.1

0.0

Difference

0.6 0.4

−0.3

−0.2

0.2

VND from nested optimization

0.8

0.3

a

technique (kinetic modeling, graphical analysis, etc.) and these estimates, together with the estimate of VND can be used to estimate any of the common outcome measures. Based on our experience, best estimates of VND can be expected when using a relatively small number of ROIs that exhibit a range of kinetic behavior. Because the signal-to-noise ratio is quite low at the voxel level, using voxel-level TACs for estimating VND cannot be expected to provide reliable estimates. An alternative to the approach described above is to regard the problem as a single (3R + 1)-dimensional optimization, i.e., to seek to optimize Eq. (6) directly. One useful general-purpose tool for highdimensional optimization problems is the simulated annealing algorithm (Kirkpatrick et al., 1983), which involves iteratively taking random steps through the optimization space, gradually reducing the range of the possible steps until convergence. This approach helps to guard against getting “stuck” in a local minimum and thus can give a greater likelihood of identifying the global minimum. In a similar context, Wong et al. (2002) proposed to use this algorithm for estimating the arterial input function CP while fitting several regions simultaneously, and Ogden et al. (2010) explored the empirical performance of this approach in several radioligands. The similarity in these procedures and the problem considered here is that in both situations, some information (either about the function CP or about the constant VND) common to the entire brain is contained within the data for each region, and that by fitting many regions at once, this information can be combined. Relative to the approach described earlier, implementing a simulated annealing algorithm is somewhat more involved, as several tuning parameters must be chosen (initial temperature, temperature reduction

0.1

0.2

0.3

0.4

0.5

VND from 1TC in the RR

0.6

0.1

0.2

0.3

0.4

0.5

0.6

VND from 1TC in the RR

Fig. 4. (Left panel) The x-axis represents the estimated VND obtained by fitting a one-tissue kinetic model to the reference region, consisting of cerebellar white matter. The y-axis represents VND estimated from the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

b 0.5 0.0 −0.5

Difference

−1.0

8 6

−1.5

4 0

−2.5

−2.0

2

VT from nested optimization

10

a

239

0

2

4

6

8

10

0

12

2

4

6

8

10

12

VT from 2TCC

VT from 2TCC

Fig. 5. (Left panel) The x-axis represents the estimated VT obtained by replacing rate constants in Eq. (1) with corresponding estimates from fitting a two tissue model to each region constraining K1/k2 to be equal to corresponding estimates from fitting a two tissue model to each region constraining K1/k2 to be equal to that estimated using a one-tissue kinetic model for the reference region, consisting of cerebellar white matter. The y-axis represents VT estimated by Eq. (1) using the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

factor, initial step sizes for each parameter, etc.), and the performance of estimation may depend on these choices. Some may prefer the nested approach because it is deterministic, whereas the simulated annealing algorithm depends on randomly determined steps. We implemented a simulated annealing algorithm for the problem of estimating VND without a reference region and repeated the analysis on both the simulated data and the real data described here (results not shown). The results are generally satisfactory, though estimation tends not to be as accurate as that for the nested optimization approach. In addition, the simulated annealing algorithm carries much greater computational demand.

b

0.0 −0.5 −1.5

−1.0

Difference

6 4 0

−2.0

2

BPP from nested optimization

8

0.5

10

a

Alternative approaches to optimization which we have not explored include Markov Chain Monte Carlo methods and particle swarm optimization (Poli et al., 2007). One choice that must be made is how many regions (and which ones) to use in the fitting. It is not necessary to include all regions, since VND can be estimated from any R TACs and then this estimate can be used to estimate binding measures for all other regions (or for all voxels). Since every region contains some information about VND, it would seem sensible to use as many as are available, but much of the information may be redundant, and more regions can increase the

0

2

4

6

8

BPP from 2TCC

10

12

0

2

4

6

8

10

12

BPP from 2TCC

Fig. 6. (Left panel) The x-axis represents the estimated BPP obtained by Eq. (5) in which VND is estimated by fitting a one-tissue kinetic model to the reference region, consisting of cerebellar white matter. The y-axis represents BPP estimated by Eq. (5) using the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

0.2 0.0

3.5

4.0

Difference

0.4

4.5

5.0

b

3.0

VND from nested optimization

a

0.6

240

3.0

3.5 4.0 VND from Lassen plot

3.0

4.5

3.5 4.0 VND from Lassen plot

4.5

Fig. 7. (Left panel) The x-axis represents the estimated VND obtained by the Lassen plot. The y-axis represents VND estimated from the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

complexity of the optimization problem and the computational expense.. With the approach described here, doubling the number of regions will require twice the computational time but will not have an effect on the stability of the optimization problem. Our experience has been that choosing several regions with a variety of kinetic behavior gives reasonably stable estimates of VND. In contrast to approaching the problem as a high-dimensional optimization problem, the nested optimization algorithm can be considered in terms of a one-dimensional optimization, albeit more complex because of the “sub-optimization” steps required. Results from direct application of standard one-dimensoinal optimization algoriothms were

b

Difference

−3

−2

−1

0

25 20 15

−5

5

−4

10

BPP from nested optimization

1

30

2

a

generally disappointing; we have found that estimation of VND tends to be more stable when we perform a simple (nested) grid search. Profile plots like those in Fig. 2 can be smooth, but for some scan datasets they can also contain local small-scale features that may induce a standard optimization routine into converging to a local minimum. As suggested by a reviewer, and confirmed by our investigations, this results from non-convexity of at least one of the ROI-specific objective functions. The relationship between VND and the estimated rate parameters for each region is generally smooth, but in some regions and some parameters there are marked discontinuities. This suggests that, over the range of VND considered, there may be identifiability issues, and we

5 10 15 20 25 30 BPP from LEGA ( VND from Lassen plot)

5 10 15 20 25 30 BPP from LEGA ( VND from Lassen plot)

Fig. 8. (Left panel) The x-axis represents the estimated VT obtained by applying LEGA. The y-axis represents VT estimated by Eq. (1) using the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

b

0 −1 −4

−3

−2

Difference

20 15 10

−5

5

BPP from nested optimization

1

25

a

241

5

10

15

20

25

BPP from LEGA ( VND from Lassen plot)

5

10

15

20

25

BPP from LEGA ( VND from Lassen plot)

Fig. 9. (Left panel) The x-axis represents the estimated BPP obtained by Eq. (5) in which VT is estimated using LEGA and VND is estimated by the Lassen plot. The y-axis represents BPP estimated by Eq. (5) using the proposed method without assuming the existence of a reference region. The solid line represents identity. (Right panel) The corresponding Bland-Altman plot.

have seen this with, e.g., k3 and k4 where both parameters are discontinuous as a function of VND but the ratio k3/k4 is relatively stable. In any case, we are interested in selecting the value of VND that gives best overall fit. Thus, we deem it worth the additional required computational expense to use a grid search that will identify a global minimum with high probability. As for computational expense, the nested optimization algorithm requires one nonlinear estimation for each ROI considered for each possible value of VND considered. For the grid search, this will be fixed by design (100 values would be a fairly fine grid), and when using the one-dimensional optimization routine, convergence generally occurs with fewer than 10 iterations. Thus with a grid of 100 values with 5 ROIs considered, an additional 500 nonlinear estimations would be required. The simulated annealing algorithm is more computationally demanding. While it depends on various settings of the algorithm, we find that for fitting one subject, the grid search implementation runs in about 5 minutes and the simulated annealing implementation runs in about 20 minutes on a modern desktop machine. As suggested by a reviewer, an alternative method for estimating a VND value that is assumed to be common across all regions is simply to fit a two-tissue model to all regions, obtain an estimate of VND (from the estimated parameters K1 and k2) for each, and then to calculate a weighted average of these. This is certainly a sensible alternative, and is likely to give reasonable estimates, at least for some radioligands. Our general approach is to minimize the total (across regions and across time points) sum of squares (equivalently, to maximize a likelihood assuming a Gaussian distribution for the errors). Thus, the choice of VND that minimizes (6) will be the one that gives best fit overall (specifically for the observed data); any other choice for VND will give larger overall mean square error for these data (and presumably, worse estimates of region-specific parameters). In some situations, the weighted mean will not be reliable—particularly cases in which parameters in a full two-tissue model are not well identified. Thus, extreme estimates of VND may result just due to identifiability issues, possibly influencing the resulting estimate of VND. This could be at least partly ameliorated by considering more robust measures of the mean, however. To study how this approach compared with those described here, we estimated VND based on a weighted average of region-specific estimates of VND obtained by fitting an unconstrained two tissue model to

the same data displayed in Figs. 4-6. Results were disappointing, as estimates had a substantial positive bias (relative to the estimate of VND based on the known reference region), and they were also considerably more variable than estimates obtained from procedures outlined here (data not shown). Thus, we tend to favor simultaneous estimation approaches for estimating a common VND value, even with the additional computational complexity. This manuscript has focused on the situation in which arterial blood data gathered during the scanning session allows for estimation of the plasma input function and thus full quantification is possible. Arterial sampling is done routinely by some groups, but it is invasive and requires significant effort to gather and model the data. In the absence of arterial blood data, there are a variety of methods available for estimating the input function (see Zanotti-Fregonara et al., 2011 for a review). The methods for estimating VND described here can be applied as long as an estimate of the input function is available, regardless of the methods used for obtaining it. Given the similarity of the simulated annealing approach for estimating VND proposed here and simulated annealing based approaches for estimating the input function (Wong et al., 2002; Ogden et al., 2010), in the absence of both a reference region and an input function, it would seem natural to estimate everything in a single expanded simulated annealing algorithm. Our recommendation in such a case, however, is to perform these two estimations in sequence, i.e., estimate the input function, and then, using the estimated input function, apply either method described here to estimate VND. Our reasoning for this is that the input function can be estimated without knowing VND, so there is no need to combine all parameters in a single optimization. Another point to make in the evaluation of the plots in Results is that, even though the reference region for the [11C]WAY-100635 radioligand has been well validated, the estimates resulting from applying a one-tissue kinetic model to the reference region do not necessarily represent “ground truth”. Although the reference region in this case is relatively large and the observations are less variable than those made on many other ROIs, there is still noise present, and so the resulting estimates are naturally affected. In addition, these estimates are based on a single region while estimates through simultaneous estimation are based on common information contained in data from several regions.

242

R. Todd Ogden et al. / NeuroImage 108 (2015) 234–242

Eliminating the requirement for a reference region will allow more radioligands to be feasible in clinical and research studies. As long as at least three of the parameters in the two-tissue model are identifiable, it will be possible to unbiasedly estimate outcome measures of interest in an automated objective manner. However, for radioligands for which only two parameters are identifiable, this leads to instability of estimation of VND using the methods described here. Such is the case for the serotonin transporter radioligand [11C]DASB (Wilson et al., 2000, 2002; Houle et al., 2000; Frankle et al., 2004) which can generally be modeled satisfactorily using a one tissue model (Ogden et al., 2007) (data not shown). Generally for this radioligand, the estimation of VND is not stable, and thus the method described here is not applicable. Disclosures/Conflicts of interest None of the authors have any disclosures to make or any conflicts of interest to report. References DeLorenzo, C., Kumar, J.S.D., Zanderigo, F., Mann, J.J., Parsey, R.V., 2009. Modeling considerations for in vivo quantification of the dopamine transporter using [11C]PE2I and positron emission tomography. J. Cereb. Blood Flow Metab. 29, 1332–1345. Forster, E.A., Cliffe, I.A., Bill, D.J., et al., 1995. A pharmacological profile of the selective silent 5-HT1A receptor antagonist, WAY-100635. Eur. J. Pharmacol. 281, 81–88. Frankle, W.G., Huang, Y., Hwang, D.R., Talbot, P.S., Slifstein, M., Van Heertum, R., Abi-Dargham, A., Laruelle, M., 2004. Comparative evaluation of serotonin transporter radioligands 11CDASB and 11C-McN 5652 in healthy humans. J. Nucl. Med. 45, 682–694. Hirvonen, J., Kajander, J., Allonen, T., Oikonen, V., Någren, K., Hietala, J., 2007. Measurement of serotonin 5-HT1A receptor binding using positron emission tomography and [carbonyl-11C]WAY-100635: considerations on the validity of cerebellum as a reference region. J. Cereb. Blood Flow Metab. 27, 185–195. Houle, S., Ginovart, N., Hussey, D., Meyer, J.H., Wilson, A.A., 2000. Imaging the serotonin transporter with positron emission tomography: initial human studies with [11C] DAPP and [11C]DASB. Eur. J. Nucl. Med. 27, 1719–1722. Huesman, R.H., Coxson, P.G., 1997. Consolidation of common parameters from multiple fits in dynamic PET data analysis. IEEE Trans. Med. Imaging 16, 675–683. Innis, R.B., Cunningham, V.J., Delforge, J., et al., 2007. Consensus nomenclature for in vivo imaging of reversibly binding radioligands. J. Cereb. Blood Flow Metab. 27, 1533–1539. Kirkpatrick, S., Gelatt Jr., C.D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671–680. Lassen, N.A., Bartenstein, P.A., Lammertsma, A.A., et al., 1995. Benzodiazepine receptor quantification in vivo in humans using [11C]flumazenil and PET: application of the steady-state principle. J. Cereb. Blood Flow Metab. 15, 152–165. Milak, M.S., Severance, A.J., Ogden, R.T., et al., 2008. Modeling considerations for 11 C-CUMI-101, an agonist radiotracer for imaging serotonin 1A receptor in vivo with PET. J. Nucl. Med. 49, 587–596.

Milak, M.S., DeLorenzo, C., Zanderigo, F., Prabhakaran, J., Dileep Kumar, J.S., Majo, V.J., Mann, J.J., 2010. In vivo quantification of human serotonin 1A receptor using 11C-CUMI-101, an agonist PET radiotracer. J. Nucl. Med. 51, 1892–1900. Milak, M., Severance, A., Prabhakaran, J., et al., 2011. In vivo serotonin-sensitive binding of [11C]CUMI-101: a serotonin 1A receptor agonist positron emission tomography radiotracer. J. Cereb. Blood Flow Metab. 31, 243–249. Millet, P., Graf, C., Moulin, M., Ibáñez, V., 2006. SPECT quantification of benzodiazepine receptor concentration using a dual-ligand approach. J. Nucl. Med. 47, 783–792. Ogden, R.T., 2003. Estimation of kinetic parameters in graphical analysis of PET imaging data. Stat. Med. 22, 3557–3568. Ogden, R.T., Ojha, A., Erlandsson, K., Oquendo, M.A., Mann, J.J., Parsey, R.V., 2007. In vivo quantification of serotonin transporters using [11C]DASB and positron emission tomography in humans: modeling considerations. J. Cereb. Blood Flow Metab. 27, 205–217. Ogden, R.T., Zanderigo, F., Choy, S., Mann, J.J., Parsey, R.V., 2010. Simultaneous estimation of input functions: an empirical study. J. Cereb. Blood Flow Metab. 30, 816–826. Parsey, R.V., Slifstein, M., Hwang, D.-R., et al., 2000. Validation and reproducibility of measurement of 5-HT1A receptor parameters With [carbonyl-llC]WAY-l00635 in humans: comparison of arterial and reference tissue input functions. J. Cereb. Blood Flow Metab. 20, 1111–1133. Parsey, R.V., Ogden, R.T., Mann, J.J., 2003. Determination of volume of distribution using likelihood estimation in graphical analysis: elimination of estimation bias. J. Cereb. Blood Flow Metab. 23, 1471–1478. Parsey, R.V., Arango, V., Olvet, D.M., Oquendo, M.A., Van Heertum, R.L., Mann, J.J., 2005. Regional heterogeneity of 5-HT 1A receptors in human cerebellum as assessed by positron emission tomography. J. Cereb. Blood Flow Metab. 25, 785–793. Parsey, R.V., Ogden, R.T., Miller, J.M., et al., 2010. Higher serotonin 1A binding in a second major depression cohort: modeling and reference region considerations. Biol. Psychiatry 68, 170–178. Poli, R., Kennedy, J., Blackwell, T., 2007. Particle swarm optimization. Swarm Intell. 1, 33–57. Raylman, R.R., Hutchins, G.D., Beanlands, R.S.B., Schwaiger, M., 1994. Modeling of carbon11-acetate kinetics by simultaneously fitting data from multiple ROIs coupled by common parameters. J. Nucl. Med. 35, 1286–1291. Turkheimer, F.E., Edison, P., Pavese, N., et al., 2007. Reference and target region modeling of [11C]-(R)-PK11195 brain studies. J. Nucl. Med. 48, 158–167. Wilson, A.A., Ginovart, N., Schmidt, M., Meyer, J.H., Threlkeld, P.G., Houle, S., 2000. Novel radiotracers for imaging the serotonin transporter by positron emission tomography: synthesis, radiosynthesis, and in vitro and ex vivo evaluation of (11)C-labeled 2(phenylthio)araalkylamines. J. Med. Chem. 43, 3103–3110. Wilson, A.A., Ginovart, N., Hussey, D., Meyer, J., Houle, S., 2002. In vitro and in vivo characterisation of [11C]-DASB: a probe for in vivo measurements of the serotonin transporter by positron emission tomography. Nucl. Med. Biol. 29, 509–515. Wong, K.-P., Meikle, S.R., Feng, D., Fulham, M.J., 2002. Estimation of input function and kinetic parameters using simulated annealing: application in a flow model. IEEE Trans. Nucl. Sci. 49, 707–713. Yaqub, M., van Berckel, B.N., Schuitemaker, A., et al., 2012. Optimization of supervised cluster analysis for extracting reference tissue input curves in (R)-[(11)C]PK11195 brain PET studies. J. Cereb. Blood Flow Metab. 32, 1600–1608. Zamuner, S., Rabiner, E.A., Fernandes, S.A., et al., 2012. A pharmacokinetic PET study of NK1 receptor occupancy. Eur. J. Nucl. Med. Mol. Imaging 39, 226–235. Zanotti-Fregonara, P., Chen, K., Liow, J.-S., Fujita, M., Innis, R.B., 2011. Image-derived input function for brain PET studies: many challenges and few opportunities. J. Cereb. Blood Flow Metab. 31, 1986–1998.