- Email: [email protected]
Likelihood estimation of drug occupancy for brain PET studies
NeuroImage 178 (2018) 255–265
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/neuroimage
Likelihood estimation of drug occupancy for brain PET studies Martin Schain a, *, Francesca Zanderigo a, b, R. Todd Ogden a, b, c a
Department of Psychiatry, Columbia University, New York, NY, USA Molecular Imaging and Neuropathology Division, New York State Psychiatric Institute, New York, NY, USA c Department of Biostatistics, Mailman School of Public Health, Columbia University, New York, NY, USA b
A B S T R A C T
Neuroimaging with PET is unique in its capability to measure in vivo the occupancy of a drug. The occupancy is typically obtained by conducting PET measurements before and after administration of the drug. For radioligands for which no reference region exists, however, the only established procedure to estimate the occupancy from these data is via linear regression analysis, forming the basis for the so-called Lassen plot. There are several reasons why simple linear regression analysis is not ideal for analyzing these data, including regression attenuation and correlated errors. Here, we propose the use of Likelihood Estimation of Occupancy (LEO) in such a situation. Similar to the Lassen plot, LEO uses the total distribution volume estimates at baseline and at block condition as input, but estimates the non-displaceable distribution volume (VND) and fractional occupancy (Δ) via direct maximum likelihood estimation (MLE). This study outlines the rationale for using MLE to estimate Δ and VND from PET data, and evaluates its performance in relation to the Lassen Plot via two separate simulation experiments. Finally, LEO and Lassen plot are applied to a PET dataset acquired with [11C]WAY-100635. LEO can exploit the covariance structure of the data to improve the accuracy and precision of the estimates of Δ and VND. Theoretically, the covariance matrix can be extracted from a test-retest dataset for the radioligand at hand. Several procedures to estimate the covariance matrix were considered as part of the simulation experiments, and the effect of the test-retest sample size was also assessed. The results are conclusive in that MLE can be used to estimate Δ and VND from PET data, avoiding the limitations associated with linear regression. The performance of LEO was, naturally, dependent on the procedure used to estimate the covariance matrix, and the test-retest sample size. Given a test-retest sample size of at least 5, but preferably 10 individuals, LEO provides higher accuracy and precision than Lassen plot in the estimation of Δ and VND. We conclude that LEO is valuable in drug occupancy studies.
Introduction Neuroimaging with Positron Emission Tomography (PET) has proven to be a valuable asset for drug development. With PET, relevant characteristics of new molecular entities can be evaluated, including level of brain exposure, target engagement (occupancy), and pharmacological activity (Gunn and Rabiner, 2017; Hargreaves et al., 2015), allowing for effective screenings already at early developmental phases. Receptor occupancy, i.e., the fraction by which a drug occupies receptors in the brain, is typically estimated from two PET measurements (one taken at baseline and the other taken after administration of the drug), using a radioligand that binds specifically to the same target. Since the drug and the PET radioligand compete for the binding sites, the difference between the two scans can be used to estimate the occupancy of the drug. Further, with well-designed studies, measurement of occupancy via PET can be used to determine the dose required for the occupancy to fall within suitable therapeutic windows (i.e., the dose range required for a clinical effect, but with limited undesirable side
effects) (Farde et al., 1988). These studies typically consist of several baseline-block scan pairs, where varying doses of the drug are administered prior to the second scan. The relationship between drug plasma concentration and the corresponding receptor occupancy can then be used to determine the drug plasma concentration required to occupy 50% of the available receptors ðKi Þ, and the maximal possible occupancy of the drug (Δmax ). The procedure with which the fractional occupancy (Δ) is estimated from PET studies, however, still deserves some attention. Binding competition between the drug and the radioligand is assumed to occur only at the specific binding sites. In addition to being specifically bound to the target of interest, however, all PET radioligands experience some degree of non-specific binding, i.e., binding to other sites than the target of interest. Also, some fraction of the radioligand is always free in tissue water. The PET system alone cannot disentangle signal originating from specifically bound radioligand from the non-displaceable (¼ free þ nonspecifically bound) radioligand. Therefore, unless the non-displaceable signal can be estimated, and in turn, the amount of specifically bound
* Corresponding author. 101 Riverside Drive, New York, NY, 10032, USA. E-mail addresses: [email protected], [email protected] (M. Schain). https://doi.org/10.1016/j.neuroimage.2018.05.017 Received 17 November 2017; Received in revised form 27 March 2018; Accepted 5 May 2018 Available online 9 May 2018 1053-8119/© 2018 Elsevier Inc. All rights reserved.
M. Schain et al.
NeuroImage 178 (2018) 255–265
radioligand can be calculated, Δ is not easily accessible. For some of the available radioligands, a brain reference region with negligible concentrations of specifically bound radioligand can be designated. For such a radioligand, Δ can easily be calculated, and they are therefore not considered in this study. For radioligands for which no reference region can be designated, estimation of Δ is not straightforward. As of today, the only established procedure for this estimation is a graphical analysis approach, in which the difference between the total distribution volumes (VT) from the two scans is plotted against that of the baseline scan, and Δ and nondisplaceable distribution volume (VND) can be obtained from linear regression. This Lassen plot was initially proposed by Lassen et al. and later revisited by Cunningham et al. (Cunningham et al., 2010; Lassen et al., 1995) There are several reasons to seek alternatives to the Lassen plot. From a theoretical standpoint, the linear regression performed in such an analysis is not strictly valid. One required assumption of regression analysis is that the independent variable is observed without measurement error; when the independent variable is measured with error, this leads to the extensively studied attenuation bias (Spearman, 1904) and has been demonstrated in PET applications (Slifstein and Laruelle, 2000). Furthermore, since baseline VT estimates are a component of both the independent and the dependent variables, the errors in the independent and dependent variables are correlated, which further violates the required assumptions in regression analysis. In this study, we present an alternative concept on how to analyze PET data acquired in blocking studies using radioligands for which no reference region exists. The general idea is to find the estimates of VND and Δ that maximize the likelihood of the observed VT values, using the same assumptions as those implicitly encoded into the Lassen plot, which are: 1) Δ is the only parameter that differs between the two scans, and 2) every region of interest (ROI) that is included in the analysis has the same value for Δ and VND. With these assumptions, joint probability density functions for the two scans can be formulated, and its maximum can be determined using standard likelihood techniques. This study is divided into 4 sections. First we provide formal mathematical derivation of the proposed algorithm, with particular emphasis on how to account for the covariance structure of the data. Second, we evaluate the algorithm's performance on simulated [11C] DASB data, including a large number of possible scenarios (varying levels of true occupancy, noise levels in the VT estimates, and procedures to estimate the covariance matrix). Third, we simulate a scenario in which estimation of Ki of a drug is sought. Finally, we test the algorithm's performance on a real dataset, in which [11C]WAY-100635 measurements were conducted before and after administration of pindolol, a partial agonist for the serotonin 1A (5HT1A) receptor. In all these scenarios, the results from the new method are compared to those obtained from the Lassen plot.
of the observed data can be expressed as
univariate Gaussian distributions, with means VS;j þ VND and pffiffiffiffiffiffi ð1 ΔÞVS;j þ VND , respectively, and standard deviation Σj;j . Further, baseline block ; VT Þ denote the joint probability density let L ðΔ; VND ; VS VT function for Vbaseline and Vblock . The analytical expression for the logT T likelihood function ℓ, neglecting the constant term that does not depend on observed data or unknown parameters, is as follows:
T
ℓ Δ; VND ; VS jVbaseline ; Vblock VND 1 VS Σ1 Vbaseline ¼ Vbaseline T T T T
VND 1 VS T VND 1 ð1 ΔÞVS Σ1 Vblock Vblock T T VND 1 ð1 ΔÞVS (3) As a function of the free parameters, equation (3) can be thought of as a surface in a kþ2 dimensional space. The coordinates at which ℓ reaches a maximum thus represent the optimal estimate for VND, Δ, and the k estimates of VS (one per ROI). Direct maximization over high dimensional space can be challenging, but in this case, it is possible to apply profile likelihood techniques and greatly reduce the dimensionality of the problem. To maximize (3) we can solve ∂V∂ S ℓ ¼ 0. After some rearrangements, we find that VS ¼
Vbaseline VND þ ð1 ΔÞ Vblock VND T T 1 þ ð1 ΔÞ2
:
(4)
For clarity purposes, we define 8 baseline > ¼ Vbaseline VND T < VS > :
Vblock ¼ S
Vblock VND T 1Δ
;
(5)
and equation (4) can be expressed as VS ¼
Vbaseline þ ð1 ΔÞ2 Vblock S S 1 þ ð1 ΔÞ2
:
(6)
Thus, given VND and Δ, there exist closed-form expressions for the maximum likelihood estimates of the other parameters, so that the loglikelihood function defined in equation (3) becomes a function of only two parameters (VND and Δ). The maximum of this function can be determined using a standard optimization routine. In our implementation, we used fminsearch. m in MATLAB R2016a. As expected, equation (6) shows that each entry of VS is replaced with the weighted average of VS estimated from the two scans, with weights determined by the occupancy. As proof of concept, Fig. 1 shows the surface spanned by equation (3) with the substitutions shown in equation b ¼ I. In addition to a global (4) in the absence of noise, and with Σ
Theory Likelihood estimation of occupancy (LEO) Let Vbaseline and Vblock denote the k 1 vectors of estimated VT values T T before and after administration of the drug in k ROIs for a given subject. According to a kinetic model with an arbitrary number of compartments, ¼ VS þ VND þ εbaseline Vbaseline T ; Vblock ¼ ð1 ΔÞVS þ VND þ εblock T
(2)
where 1 denotes a k 1 vector of ones, and Σ represents the covariance matrix of the data, whose estimation is discussed in a separate section (for now, we will treat Σ as if it is known). This means that, for ROI j, baseline block the marginal probability density functions for VT;j and VT;j are
Material and methods
Vbaseline NðVND 1 þ VS ; ΣÞ T ; Vblock NðVND 1 þ ð1 ΔÞVS ; ΣÞ T
minimum at the correct coordinates ðΔ ¼ 0:5; VND ¼ 2Þ; a curved ridge with a local minimum is observed in the parameter space, possibly reflecting that other combinations of VND and Δ may provide reasonable fits in some instances.
(1)
where Δ denotes the fractional occupancy of the drug, 0 Δ 1, VS is the distribution volume of specifically bound radioligand, and εbaseline and εblock are error terms associated with each measurement. In this model, equivalently to Lassen plot, VS varies across the k ROIs, whereas Δ and VND are constants. Assuming normality of the errors, the distribution
Estimation of covariance matrix The covariance matrix Σ in equations (2) and (3) represents the within-subject covariance matrix of the regional VT values across ROIs. This can be interpreted as the covariances across regional VT values for 256
M. Schain et al.
NeuroImage 178 (2018) 255–265
the sampled eigenvalues, and derive a procedure to non-linearly shrink the sample covariance matrix towards the population asymptote. The shrunken matrix may therefore represent a better approximation of the real, but unknown, covariance matrix than the sample covariance matrix. The procedure introduces no strong assumptions regarding the distribution of the population eigenvalues. Note that with any of the methods to estimate Σ among those outlined above (Identity, Diagonal and Shrink), the number of ROIs included in the analysis is not limited by the number of test-retest subjects. The consideration of the covariance between measurement errors in different ROIs marks a distinct difference to an earlier study that pursued a similar idea (Radua et al., 2011). In that study, the maximum likelihood formulation considered independent Gaussian noise with a fixed variance in all ROIs. In the current context, such procedure corresponds to using the identity matrix multiplied by a free parameter (σ ), whose value is estimated by the maximum likelihood estimator. However, this constant term does not affect the coordinates of the minimum of the log-likelihood function, and therefore, the previously proposed procedb ¼ I. The possibility to estimate Σ ure is in effect the same as setting Σ
Fig. 1. LEO proof-of-concept. The figure shows the surface spanned by ℓ over a range of values for Δ and VND. The figure shows ℓ rather than ℓ for visualization purposes. A minimum is observed at the correct coordinates ðΔ ¼ 0:5; VND ¼ 2Þ.
directly from test-retest data, without requiring strong assumptions of measurement noise characteristics, marks an additional distinction from the approach outlined in a recent study, which proposed to use a noise propagation model to estimate the elements of Σ (Naganawa et al., 2017).
one subject who has undergone repeated PET measurements, so that the differences in VT across scans are in large attributable to measurement errors, rather than biological variability. Since Σ contains covariances between all ROIs, it is a k k matrix of k unknown parameters (k þ parameters in the most general case). It 2 is therefore not generally feasible with two PET scans to estimate the elements of Σ in the same analysis in which we are estimating VND and Δ. Fortunately, for many radioligands, data from a test-retest experiment, in which each of several subjects are imaged twice, are available. In a testretest experiment biological variability between the test and retest scan is normally assumed to be negligible. It can be shown, that if Vtest and T Vretest denote the k N matrices, where each of the N columns contains k T regional VT values for one of the N subject included in the test-retest experiment, then the sought covariance matrix becomes Σ covðVtest T Vretest Þ=2 if N is sufficiently large (see the supplementary material for T derivation). However, unless the number of subjects included in the test-retest study is large, the covariances between ROIs may not be accurately estimated. In addition, if the number of test-retest subjects is smaller than the number of ROIs included in the analysis (i.e., N < k), the estimate of Σ will not be invertible. Most test-retest studies include 5–10 research subjects, and in such a case, we cannot count on the raw estimate of Σ being sufficiently accurate. In this study, we evaluate 3 different procedures to approximate Σ from a small-sampled test-retest study, each b with increasing complexity (estimates of Σ are denoted Σ).
Simulating data with known covariance structure It is likely that the estimation of the covariance matrix will impact the performance of LEO. To provide a “gold standard”, unattainable in practice, we will consider the situation in simulations in which the “true” covariance structure of the data is known. Let Σ denote the true, known covariance matrix for a given data set. The Cholesky factorization of Σ produces an upper triangular matrix Γ, so that Σ ¼ ΓT Γ. Let B denote a k M matrix, where each column holds k independent random numbers, sampled from Nð0; 1Þ. Then, covðΓT BÞ ¼ ΓT covðBÞΓ ¼ ΓT IΓ ¼ Σ: Let Vtrue denote a k 1 vector representing true T (i.e., noise free) VT values in k ROIs. Then, by adding the columns in ΓT B true are created, each with unique noise, and to Vtrue T , M realizations of VT with covariance matrix equal to Σ. Depending on how Vtrue is defined, T this procedure allows for simulation of both blocking and test-retest data.
Simulation of [11C]DASB blocking data In this simulation, we seek to investigate the overall performance of the various options for LEO in comparison to the Lassen plot, in terms of accuracy and precision in estimates of Δ and VND. We evaluate all 3 b (identity, diagonal, and shrink), various levels of different versions of Σ noise, a range of values for the true occupancy, and different number of b is subjects included in the hypothetical test-retest study from which the Σ obtained. A flowchart describing the simulation setup is reported in Fig. 2. This simulation is based on parameters estimated from a previously published test-retest dataset with the serotonin transporter radioligand [11C]DASB in 10 subjects (Ogden et al., 2007) (test and retest scans were was defined as the conducted on the same day). From this dataset, Vtrue T average VT across individuals and scans in k ¼ 8 ROIs that, based on previous publications, encompass a range of VT-values (Parsey et al., 2006) (anterior cingulate cortex, amygdala, dorsal putamen, insula, occipital cortex, orbitofrontal cortex, thalamus, and ventral striatum). Regional VT values were calculated using the 1-tissue compartment model, as previously suggested (Ogden et al., 2007; Ginovart et al., 2001). Further, an underlying “true” covariance matrix, Σ, was defined Vretest Þ=2, where Vtest and Vretest denote 8 10 via Σ ¼ covðVtest T T T T matrices holding regional VT values from the test and retest scans, respectively. We take this estimate of Σ to be the “true” Σ for the purposes of our simulation.
b ¼ I. Setting Σ b to equal the identity matrix corresponds to 1) Identity, Σ disregarding all differences in variance among regions and any covariance between regions. b ¼ diagðcovðVtest Vretest Þ=2Þ. Setting Σ to equal the di2) Diagonal, Σ T
T
agonal version of the estimated covariance matrix corresponds to considering distinct variances for each ROI, but setting to zero the covariances between ROIs. b obtained from non-linear shrinkage of the sample covari3) Shrink, Σ ance matrix towards the population asymptotes. Option 3 above was obtained from the work of Ledoit and Wolf, who derived a shrinkage formula for covariance matrices estimated with small sample sizes (Ledoit and Wolf, 2015). In brief, Ledoit and Wolf present a method that consistently approximates the population eigenvalues from
257
M. Schain et al.
NeuroImage 178 (2018) 255–265
above procedure. In each simulation set, 1000 instances were used (i.e., in (7), i ¼ 1; …; 1000). Parameters that varied among the different simulations were: occupancy (Δtrue ¼ 0.25, 0.50, and 0.75), noise levels (λ ¼ 1, 1.1, and 1.25, corresponding to 100%, 110% and 125% of realistic b (ranging levels), and number of test-retest subjects used to calculate Σ between 5 and 100). For each noise instance, LEO and Lassen plot were b and V d used to estimate Δ and VND (estimates denoted Δ ND ). LEO was employed with three different estimates of Σ described in 2.1.2. As a proof of concept, LEO was also employed with the Oracle covariance matrix (i.e., the matrix obtained from the measured test-retest data, which represents the true Σ in this simulation). For each simulation, and b the squared distance d between the estimated and for each choice of Σ, the dΔ ¼
true value for occupancy and VND was P1000 true P 2 b 2 ; dVND ¼ 1000 ðV true V d ΔÞ ND Þ . ND i¼1 ðΔ i¼1
calculated,
Simulation of a blocking experiment to determine Ki In this section, we simulate a scenario in which the inhibition constant Ki for a drug is sought. In this scenario, it is common that following the baseline scan, research subjects are administered different doses of the drug for the blocking scan. The baseline and blocking scans are used to estimate the Δ at each dose. With Δ estimated at different dose levels, Ki and the maximal occupancy, Δmax , can be fitted in
Fig. 2. Flow-chart describing the setup for the [11C]DASB simulation.
With Σ and Vtrue defined, the procedure outlined in the previous T section (Simulating data with known covariance structure) was used to simulate new test-retest subjects (i.e., the procedure was performed twice with the same values for Vtrue T ). These simulated test-retest subb following the procedures described in jects were then used to obtain Σ
Δ ¼ Δmax
administered before the second scan (Parsey et al., 2006). Vtrue;baseline T and Vtrue;block was defined as the average VT value in each ROI across T true and subjects. Subsequently, a Lassen plot was performed to define VND . As before, Σ, defined from the measured test-retest dataset Vtrue S described above, was used to create M noise instances with known covariance. Separate noise instances were generated for the baseline and true , Vtrue blocking data sets, ðΓT BÞbaseline and ðΓT BÞblock . Finally, with VND S , and Σ defined, the i'th realization of k simulated VT values at baseline and after block was calculated as baseline
true Vsim;baseline ¼ Vtrue þ VND 1 þ ðΓT BÞi S T;i
Vsim;block T;i
true
¼ ð1 Δ
Þ
Vtrue S
þ
true VND 1
λ
block
þ ðΓT BÞi
λ
;
(8)
where Cp is the concentration of the drug in blood plasma at the time of each of the block scans. Here, we define a hypothetical drug that binds specifically to the 5HT1A receptor, has a Ki of 25 ng/mL, and potential to occupy all receptors (Δmax ¼ 1). In our design, a total of 10 blocking experiments are conducted at varying doses of the drug. At the time of blocking scans, Cp was set to 5, 10, 20, 30, 40, 50, 60, 70, 80, and 90 ng/mL. Insertion of these numbers into equation (8) provides the “true” values for Δ that corresponds to the defined Cp. The procedure to simulate the PET data was similar to that outlined in the previous section. First, VT-values from a previous same-day testretest study (Parsey et al., 2000) in 5 individuals with the 5HT1A radioligand [11C]WAY-100635 were used to define Σ, and subsequently, to simulate a test-retest dataset. The true VT values were defined as the average VT across subjects and scans, in 4 arbitrarily selected ROIs associated with a range of values (anterior and posterior cingulate cortex, hippocampus, and insula). Note that, in the current analysis, the included number of ROIs had to be < 5, since only 5 test-retest subjects were available (with more ROIs, the “true” Σ would not have been invertible, which was required for noise generation). The true VND was set to 0.3, which is in the range of VND reported for [11C] WAY-10063515, and the noise scale factor λ was set to 1, corresponding to 100% of realistic noise. The procedure for this simulation, shown in Fig. 3, is as follows. First, 10 test-retest-subjects were simulated using Σ. Second, 10 baselineblocking scan pairs were simulated using equation (7), each with unique noise (also from Σ), and with Δ defined from equation (8) using b was calculated from the simulated testthe hypothetical Cp. Third, Σ
section Estimation of covariance matrix. Theoretically, the estimate of Σ will improve with increasing number of test-retest subjects. Therefore, the experiment was repeated with varying number of test-retest subjects, ranging from 5 to 30 in steps of 5, and then from 30 to 100 in steps of 10. Most brain PET test-retest studies include between 5 and 10 subjects, and therefore these settings provide information on the performance of the algorithm under realistic conditions. With larger testretest sample sizes, the errors introduced by a poorly estimated covariance matrix should hypothetically decrease, and the results gradually approach those obtained with a perfectly estimated (i.e., a true) covariance matrix. For this purpose, we also evaluated LEO's performance when the true covariance matrix was used, which is referred to as the “Oracle” method. This scenario represents an upper limit of LEO's performance, since it is insensitive to potential errors in the estimate of Σ. As such, the Oracle method enables us to distinguish between the effects of a poorly estimated covariance matrix and the uncertainty intrinsic to LEO. For the blocking data, we based our simulation on 8 previously published blocking experiments with [11C]DASB, where sertraline was
(
Cp Cp þ Ki
retest measurements, using only Shrink (the other procedures to estib were omitted in results since previous simulation was conclusive mate Σ b Fourth, VND and in that Shrink was the preferred procedure to obtain Σ). Δ were estimated both using LEO and Lassen plot. Finally, using the known Cp and the estimated values for Δ, Ki and Δmax were estimated using equation (8). Repeating this procedure 1000 times results in 1000 estimates of Ki and Δmax obtained with Lassen plot, and, 1000 estimates for LEO. The distributions of Ki and Δmax were visualized with histograms, and the agreement were quantified by calculating the squared difference
(7)
where ðΓT BÞi denote the i'th column in the noise matrix, Δtrue represents the true occupancy, and λ is a scale factor for the noise. A total of 27 sets of simulations were performed, each following the 258
M. Schain et al.
d dKi ¼
between the estimated and the true P1000 true P max 2 bi Þ2 ; dΔmax ¼ 1000 ðΔmax;true Δd ðK K Þ . i i¼1 i¼1
NeuroImage 178 (2018) 255–265
Δ was also calculated using the Lassen plot and LEO, and the results were compared to those estimated using the reference region (average across ROIs). In contrast to the simulations above, the true b was esticovariance matrix Σ is not known. When applying LEO, Σ
values,
mated from the [11C]WAY-100635 test-retest dataset in 5 individuals that was used in section Simulation of a blocking experiment to determine b was only estimated with the Shrink method, since preKi14. Here, Σ
Application to real data In this final section, we apply Lassen plot and LEO to a previously published dataset (Martinez et al., 2001), consisting of [11C] WAY-100635 measurements at baseline, and after administration of pindolol, a partial 5HT1A receptor agonist. In this dataset, arterial input functions were acquired, allowing for estimation of regional VT. In addition, for [11C]WAY-100635, the cerebellar white matter has been shown to be a suitable reference region, allowing for estimation of VND. To estimate VND, a 1-tissue compartment model was fitted to the time-activity curve for the cerebellar white matter. To estimate regional VT values, a 2-tissue compartment model was applied, with K1/k2 constrained to the estimate of VND, as previously suggested (Parsey et al., 2000). Estimates of VT and VND enables estimation of regional binding potential (Parsey et al., 2005). Here, we used BPP (¼ VT-VND) as outcome measure (Parsey et al., 2000). Using the reference region, the occupancy of pindolol was estimated via Δ ¼ ðBPbaseline P Þ=BPbaseline in 7 ROIs (amygdala, posterior cingulate cortex, BPblock P P hippocampus, occipital cortex, orbitofrontal cortex, parietal cortex, and parahippocampal gyrus). These ROIs were selected since they, after preliminary analysis of BPP, displayed similar occupancy within each subject (<10% difference from across-ROI mean).
vious analyses showed that this method is preferable over the simpler procedures to estimate Σ. Occupancies obtained using the reference region were compared to those obtained from Lassen plot and LEO, via regression analysis.
Implementation notes In the presence of noise, both LEO and the Lassen plot can “fail”, i.e., provide highly unrealistic estimates. To avoid this, constraints on the free parameters were enforced in LEO. By definition, 0 Δ 1, so the occupancy was constrained to be in this interval by performing an uncon Δ strained minimization over log 1Δ rather than over Δ. With a similar reasoning, a non-negativity constraint was enforced for VND, by minimizing over logðVND Þ. Although the Lassen plot has an implicit upper bound on the occupancy, it is generally executed without explicit constraints on the free parameters. Thus, the mere presence of parameter constraints in LEO has the potential to result in lower scores of dΔ and dVND , without actual improved accuracy. To avoid this, all occupancy and VND values obtained from the Lassen plot that were <0 were then set to 0 before dΔ and dVND were calculated. Also, the “failure rate” was defined as the number of simulation instances for which a method provided an estimate of VND that was larger than the maximal VT for that simulated subject. All simulations that failed according to this definition were excluded before calculation of dΔ and dVND . Results Simulation of [11C]DASB blocking data Results from simulation 1 are shown in Figs. 4–7. Figs. 4 and 5 report the sum of squared distances d between true and estimated values for Δ and VND for various test-retest sample sizes (thus, a low value reflects a good performance). LEO using the true underlying covariance matrix (i.e., the “Oracle” method) always resulted in the lowest squared distances. Although this is interesting from a theoretical point of view, the “Oracle” method is not possible in practice, as it requires knowledge of the true underlying covariance structure of the data. Of interest, with increasing test-retest sample sizes, using non-linear shrinkage to estimate the covariance matrix provided results that asymptotically approached those obtained with the Oracle method. Considering commonly used test-retest sample sizes (i.e., 10 or 5 subjects), the results indicate that LEO performed with covariance matrices estimated with either the Shrink method, or diagonal matrices, is preferable over Lassen plot, as these methods consistently resulted in lower scores of dΔ and dVND . The difference between LEO and Lassen plot b was less evident when the identity matrix was used in place of Σ. To allow further visualization of some of the results shown in Figs. 4 and 5, the distributions of the estimates of Δ and VND obtained with 10 test-retest subjects are reported in Figs. 6 and 7, respectively. In these violin plots, the relative width and color within each violin represent the density of the respective distributions. The black dashed line corresponds to the true value in each plot. The peaks around 0 seen for Lassen plot and LEO executed with the identity and diagonal matrix for 25% occupancy in Fig. 6 results from the lower bound constraint; all instances for which the methods would provide a negative occupancy
Fig. 3. Procedure for the simulation of a blocking experiment to determine Ki.
259
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 4. Sum of squared error in estimates of occupancy (dΔ ) for various noise levels, occupancy levels, and test-retest sample sizes used to estimate the covariance matrix Σ. The lowest error was consistently observed with the Shrink method, which asymptotically approached results obtained with the true Σ (i.e., Oracle). Shrink also resulted in smaller error than Lassen plot even for small sample sizes.
are floored at 0 (see Implementation notes). Although no clear bias was observed in occupancy estimated with the Lassen plot, there was a considerable spread. This spread was reduced with LEO when executed with Shrink approach to estimate Σ. In contrast to what was observed for occupancy, all methods resulted in a systematic overestimation of VND (Fig. 7). The largest bias was observed at 25% occupancy, with Lassen plot and LEO executed with the identity matrix in place of the covariance matrix. For higher occupancy levels, the bias was reduced. In all cases, LEO executed with the Shrink approach to approximate Σ resulted in lower bias and lower spread than the other procedures available in reality. In supplementary figure S1, the true and approximated covariance matrices are visualized for varying number of test-retest subjects. Inspection of these matrices reveals that Σ is not approximately diagonal, but rather, many off-diagonal elements are not negligible. Further, figure S2 in the supplementary material shows associations of the diagonal elements of Σ with average and standard deviation of the VT values (calculated from the baseline data), and with ROI volumes. These associations were r2 ¼ 0.98 for average VT values, r2 ¼ 0.93 for their standard deviations, and r2 ¼ 0.60 for ROI volumes (Figure S2, supplementary material).
Lastly, Table S1 in the supplementary material reports the rate by which each of the evaluated methods “fails” according to the criterion described in Implementation notes. For all methods, highest failure rates were observed at low occupancy levels (for instance, at 100% noise, Lassen plot failed in 5% of the instances at 25% occupancy, and never at 75% occupancy). Shrink, regardless of test-retest sample size, consistently resulted in the lowest failure rate (with 10 test-retest subjects and 100% noise, the failure rate was 1.4% at 25% occupancy). The highest b failure rate was observed when the identity matrix was used in place of Σ, whereas Lassen plot and LEO using a Diagonal matrix resulted in similar failure rates.
Simulation of [11C]WAY-100635 blocking data to estimate Ki From the previous experiment, it was evident that using the Shrink approach to estimate the covariance matrix in LEO provided better results than when the diagonal or identity matrix was used. For this purpose, these simplified versions were omitted from the remaining experiments. The “Oracle” method was also excluded since it is not available in real data applications.
260
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 5. Sum of squared error in estimates of VND (/1000) (dVND ) for various noise levels, occupancy levels, and test-retest sample sizes used to estimate the covariance matrix Σ. The lowest error was consistently observed in correspondence of the Shrink method, which asymptotically approached results obtained with the true Σ (i.e., Oracle). Shrink also resulted in smaller error than Lassen plot even for small sample sizes.
were r2 ¼ 0.99, slope ¼ 1.02, and intercept ¼ 0.04. The Lassen plots carried out to estimate the occupancies shown in Fig. 9 are presented in Figure S3 in the supplementary material.
Following the procedure outlined in Fig. 3 resulted in 1000 datasets, each with 10 estimated values for Δ corresponding to the defined drug plasma levels. The estimated Δ were then used in equation (8) to estimate Ki and Δmax . With Ki and Δmax estimated, the occupancy plot describing the relationship between drug plasma levels and receptor occupancy were obtained. The top panel in Fig. 8 display these relationships with occupancy estimated with LEO and the Lassen plot, respectively, and the bottom panels show the distributions of Ki and Δmax . LEO resulted in a lower spread in the estimates of Ki and Δmax compared with the Lassen plot. The squared distances for the two parameters were dKi ¼ 211:5; and dΔmax ¼ 0:034 for Lassen plot, and dKi ¼ 59:0 and dΔmax ¼ 0:010 for LEO.
Discussion Graphical procedures for the analysis of PET data are attractive, due to their simplicity to implement and interpret (Logan et al., 1990; Patlak et al., 1983). The Lassen plot is one of the most commonly used approaches to estimate occupancy using radioligands for which no suitable reference region exist, in spite of it being not strictly valid due to correlated errors and errors in the predictor variable. To avoid the assumptions implicit in linear regression, we here consider maximum-likelihood (ML). ML is a well-suited method for the problem at hand, as it relies on a solid theoretical framework assuring that the estimates of the free parameters of interest are approximately unbiased. In addition, the ML estimator is consistent, meaning that it is capable of determining the true parameter value with arbitrarily high precision given a sufficient number of measurements. In the derivation of a ML estimator for occupancy and VND, it is important not to introduce additional assumptions regarding the underlying theory of the pharmacokinetic properties of the radioligand or the drug. For this purpose, we formulated the likelihood function based
Application to real data As a final experiment, we applied LEO and Lassen plot to a previously published dataset, in which 9 subjects had [11C]WAY-100635 PET measurements obtained before and after administration of pindolol. The occupancy calculated using LEO (with Shrink) or the Lassen Plot was compared to that obtained from using BPP. The scatter plot comparing these occupancy levels is shown in Fig. 9. The linear regression analysis resulted for Lassen plot in r2 ¼ 0.83, slope ¼ 1.33, and intercept ¼ 0.20. For LEO, the corresponding values 261
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 6. Violin plots showing the distribution of occupancy estimated with Lassen plot and LEO, assuming 10 test-retest subjects. The relative width and color in each violin represents the density, and the dashed black line corresponds to the true value. No evident bias was observed for any method, but LEO executed with Shrink resulted in the lowest spread. The additional yellow fields seen at 25% occupancy result from non-negativity constraints (see section Implementation notes).
Although a perfectly acquired test-retest dataset may provide a theoretical possibility to calculate the covariance matrix, there are practical limitations. The accuracy of the estimated covariance matrix depends on the sample size of the test-retest dataset. From a test-retest sample size of 5–10 individuals, simply calculating Σ using the definition for covariance is not recommendable, as the covariances will be poorly estimated. In addition, if the number of ROIs included in the analysis is larger than the number of test-retest subjects, Σ will not be invertible, which disables application of LEO. Rather, alternative methods to estimate Σ need to be undertaken. Note that, with such methods to estimate Σ (for instance, using the diagonal matrix or nonlinear shrinkage), the number of ROIs included in the analysis is not constrained by the test-retest sample size. Here, we evaluated three such methods, and the results show that estimation of the covariance matrix is of importance. First, if equal variances and zero covariances were assumed (i.e., the identity matrix was used in place of Σ), LEO resulted in similar spread and bias as that observed for the Lassen plot. Second, if different within ROI variances were considered but covariances between ROIs were still neglected (i.e., using the diagonal matrix), a somewhat improved estimation of Δ and VND was observed. The performance improved even further when the full covariance matrix was estimated via non-linear shrinkage. Finally, in the hypothetical situation that the true covariance matrix was completely known, VND and occupancy were estimated with high accuracy and precision. This latter option can be thought of as the upper bound of the LEO's performance, which is very unlikely to ever be reached in a real application. For real clinical investigations, the results from this study indicate that the non-linear shrink method to estimate the covariance matrix based on 10 test-retest subjects is favorable in comparison the Lassen plot.
on the standard assumptions in compartmental modeling: each molecule of the radioligand that has penetrated the brain is, at any given time point, either specifically bound to the target, or non-displaceable (Innis et al., 2007). Further, the error term was assumed to follow a normal distribution. This latter assumption is in large supported by the fact that LEO takes as input the estimated VT-values, which are typically derived from non-linear regression. According to likelihood theory, it then follows that the errors in VT are (at least asymptotically) normally distributed. However, the error associated with an estimate of VT may come from multiple sources, for instance statistical noise in the emission and transmission data, error in the image reconstruction, and measurement error of the arterial input function. All of these error sources have the possibility to influence different brain regions similarly. It is therefore not reasonable to assume that the errors in two separate brain regions from the same subject are independent. To allow for dependency in errors across brain regions, the covariance structure across brain regions was considered in the derivation of LEO. In a blocking experiment, each subjects VT values are compared at baseline and after administration of the drug, and the fractional occupancy, which is treated as a free parameter, is the only parameter that is assumed to differ between the scans (i.e., VS and VND are assumed equal between baseline and block). The covariance matrix should therefore only contain covariance between measurement errors, without the influence of biological variability that is present across subjects. We here propose to approximate this information using test-retest datasets for the radioligand, since the difference between a test and a retest scan in one individual should in large reflect only measurement errors. However, under certain circumstances, some biological variability between test and retest can be expected due to, for instance, diurnal variation (Matheson et al., 2015). 262
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 7. Violin plots showing the distribution of VND estimated with Lassen plot and LEO, assuming 10 test-retest subjects. The relative width and color in each violin represents the density, and the dashed black line corresponds to the true value. For 25% occupancy, an evident overestimation in VND was observed for all methods (except the Oracle method), and the bias was reduced at higher occupancy levels. Of the methods that are typically available in real studies, LEO executed with Shrink provided the smallest bias and spread.
at hand. That LEO explicitly uses test-retest information possibly enables estimation of Δ at lower values. Last, we investigated the agreement between estimates of occupancy obtained using a reference region to those obtained using Lassen plot and LEO. As seen in Fig. 9, occupancy estimated with LEO was on group level in closer agreement with that estimated using a reference region than estimates obtained with Lassen plot. These results should however be treated with some caution, since the sample is rather small, and it is possible that, at least in part, the improved performance seen with LEO is driven by the parameter constraints. When comparing LEO and Lassen plot using real data, no particular attention was given to estimates of VND. For [11C]WAY-100635, VND is very low, typically around 0.3 when estimated using the VT in the cerebellar white matter (Parsey et al., 2005). For 5 out of the 9 subjects, Lassen plot resulted in negative values of VND (figure S3 in supplementary material), whereas with LEO these values were positive (possibly as a result from the non-negativity constraint), but very close to 0. From these observations, we considered [11C]WAY-100635 VND estimated with either method to be biologically not meaningful. Regardless of whether Lassen plot or LEO is used, it is important to be aware that both algorithms assume each ROI included in the analysis to have the same value for VND and occupancy. It is not obvious how to verify this assumption using observed data. We note that a violation of this assumption can result in an unpredictable bias that may confound findings. In the Lassen plot, a non-linear relationship or clear outliers may give hints as to whether certain ROIs should be excluded from the analysis. We suggest that when estimating parameters using LEO, that also visually examining the Lassen plot can
An interesting question is how well the covariance matrix can be estimated in a scenario where test-retest data are not available. Supplementary figure S2 shows that, using only single scan data, the diagonal elements of Σ can be very well predicted using VT values and ROI volumes, indicating that single scan data can likely be used to estimate within-ROI variances. However, as seen in supplementary figure S1, Σ is not approximately diagonal – many off-diagonal elements have nonnegligible values. Further research is needed to understand how to better approximate off-diagonal elements of Σ either using single scan data, or, as recently proposed, modeling the noise characteristics of the PET data (Naganawa et al., 2017). To assess the impact of using LEO in real life applications, we simulated a scenario in which the properties of an unknown new drug is sought, with the underlying hypothesis that the estimation Ki and Δmax will improve as a result of improved accuracy and precision of Δ. From the results in Simulation 1, the largest difference between Lassen plot and LEO was observed at low or intermediate occupancy values, where LEO provided higher accuracy and precision in the parameter estimates. Poorly estimated occupancies at low drug plasma levels (typically corresponding to low occupancies) impacts the entire fit of equation (8) and may thus produce uncertainties in the estimates of Ki, regardless of the quality of the estimates at higher plasma levels. We therefore hypothesized that, with LEO, the estimation of Ki will improve as a result of improved accuracy and precision of Δ, in particular at low drug plasma levels. The results presented in Fig. 8 confirmed this hypothesis, where a lower spread across noise instances is observed in the occupancy curves (top panel). In general, it is likely that the lower limit of a detectable occupancy is determined by the test-retest variability of the radioligand
263
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 8. Results from the simulated [11C]WAY-100635 experiment. Top panel show the relationship between receptor occupancy and plasma concentration obtained with LEO and Lassen plot (each line represents one simulation). Bottom panels show the distributions of estimates of Ki and Δmax . In each plot, the black line corresponds to the true value.
Funding This study was supported by the Paul Janssen Fellowship in Translational Neuroscience. Acknowledgements We thank Prof. Michael Wolf and Dr. Olivier Ledoit at the Department of Economics, University of Zurich, Switzerland, for helpful discussions regarding estimation of covariance matrices from small datasets. Appendix A. Supplementary data Supplementary data related to this article can be found at https://doi. org/10.1016/j.neuroimage.2018.05.017. References Cunningham, V.J., Rabiner, E.A., Slifstein, M., Laruelle, M., Gunn, R.N., 2010. Measuring drug occupancy in the absence of a reference region: the Lassen plot re-visited. J. Cereb. Blood Flow. Metab. 30, 46–50. Farde, L., Wiesel, F.A., Halldin, C., Sedvall, G., 1988. Central D2-dopamine receptor occupancy in schizophrenic patients treated with antipsychotic drugs. Arch. Gen. Psychiatry 45, 71–76. Ginovart, N., Wilson, A.A., Meyer, J.H., Hussey, D., Houle, S., 2001. Positron emission tomography quantification of [(11)C]-DASB binding to the human serotonin transporter: modeling strategies. J. Cereb. Blood Flow. Metab. 21, 1342–1353. Gunn, R.N., Rabiner, E.A., 2017. Imaging in central nervous system drug discovery. Semin. Nucl. Med. 47, 89–98. Hargreaves, R.J., Hoppin, J., Sevigny, J., et al., 2015. Optimizing central nervous system drug development using molecular imaging. Clin. Pharmacol. Ther. 98, 47–60. 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. 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.
Fig. 9. Scatter plot between occupancies obtained using a reference region to calculate BPp, and those obtained from VT values using Lassen plot and LEO. Σ was estimated using the Shrink method applied to VT values from 5 testretest subjects.
provide a good complementary analysis and potentially reveal problems with the data. Finally, to facilitate the application of LEO for other researchers; we freely distribute MATLAB code for doing the calculations. The code and instructions are available at https://github.com/martinschain/LEO. 264
M. Schain et al.
NeuroImage 178 (2018) 255–265 humans: comparison of arterial and reference tisssue input functions. J. Cereb. Blood Flow. Metab. 20, 1111–1133. Parsey, R.V., Arango, V., Olvet, D.M., Oquendo, M.A., Van Heertum, R.L., John Mann, J., 2005. Regional heterogeneity of 5-HT1A receptors in human cerebellum as assessed by positron emission tomography. J. Cereb. Blood Flow. Metab. 25, 785–793. Parsey, R.V., Kent, J.M., Oquendo, M.A., et al., 2006. Acute occupancy of brain serotonin transporter by sertraline as measured by [11C]DASB and positron emission tomography. Biol. Psychiatry 59, 821–828. Patlak, C.S., Blasberg, R.G., Fenstermacher, J.D., 1983. Graphical evaluation of blood-tobrain transfer constants from multiple-time uptake data. J. Cereb. Blood Flow. Metab. 3, 1–7. Radua, J., Bullich, S., Lopez, N., Catafau, A.M., 2011. Restricted maximum likelihood estimation of PET neuroreceptor occupancy in the absence of a reference region. Med. Phys. 38, 2558–2562. Slifstein, M., Laruelle, M., 2000. Effects of statistical noise on graphic analysis of PET neuroreceptor studies. J. Nucl. Med. 41, 2083–2088. Spearman, C., 1904. The proof and measurement of association between two things. Am. J. Psychol. 15, 72–101.
Ledoit, O., Wolf, M., 2015. Spectrum estimation: a unified framework for covariance matrix estimation and PCA in large dimensions. J. Multivar. Anal. 139, 360–384. Logan, J., Fowler, J.S., Volkow, N.D., et al., 1990. Graphical analysis of reversible radioligand binding from time-activity measurements applied to [N-11C-methyl](-)-cocaine PET studies in human subjects. J. Cereb. Blood Flow. Metab. 10, 740–747. Martinez, D., Hwang, D., Mawlawi, O., et al., 2001. Differential occupancy of somatodendritic and postsynaptic 5HT(1A) receptors by pindolol: a dose-occupancy study with [11C]WAY 100635 and positron emission tomography in humans. Neuropsychopharmacology 24, 209–229. Matheson, G.J., Schain, M., Almeida, R., et al., 2015. Diurnal and seasonal variation of the brain serotonin system in healthy male subjects. Neuroimage 112, 225–231. Naganawa, M., Gallezot, J.D., Rossano, S., Carson, R.E., 2017. Quantitative PET imaging in drug development: estimation of target occupancy. Bull. Math. Biol. Ogden, R.T., Ojha, A., Erlandsson, K., Oquendo, M.A., Mann, J.J., Parsey, R.V., 2007. In vivo quantification of serotonin transporters using [(11)C]DASB and positron emission tomography in humans: modeling considerations. J. Cereb. Blood Flow. Metab. 27, 205–217. Parsey, R.V., Slifstein, M., Hwang, D.R., et al., 2000. Validation and reproducibility of measurement of 5-HT1A receptor parameters with [carbonyl-11C]WAY-100635 in
265
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/neuroimage
Likelihood estimation of drug occupancy for brain PET studies Martin Schain a, *, Francesca Zanderigo a, b, R. Todd Ogden a, b, c a
Department of Psychiatry, Columbia University, New York, NY, USA Molecular Imaging and Neuropathology Division, New York State Psychiatric Institute, New York, NY, USA c Department of Biostatistics, Mailman School of Public Health, Columbia University, New York, NY, USA b
A B S T R A C T
Neuroimaging with PET is unique in its capability to measure in vivo the occupancy of a drug. The occupancy is typically obtained by conducting PET measurements before and after administration of the drug. For radioligands for which no reference region exists, however, the only established procedure to estimate the occupancy from these data is via linear regression analysis, forming the basis for the so-called Lassen plot. There are several reasons why simple linear regression analysis is not ideal for analyzing these data, including regression attenuation and correlated errors. Here, we propose the use of Likelihood Estimation of Occupancy (LEO) in such a situation. Similar to the Lassen plot, LEO uses the total distribution volume estimates at baseline and at block condition as input, but estimates the non-displaceable distribution volume (VND) and fractional occupancy (Δ) via direct maximum likelihood estimation (MLE). This study outlines the rationale for using MLE to estimate Δ and VND from PET data, and evaluates its performance in relation to the Lassen Plot via two separate simulation experiments. Finally, LEO and Lassen plot are applied to a PET dataset acquired with [11C]WAY-100635. LEO can exploit the covariance structure of the data to improve the accuracy and precision of the estimates of Δ and VND. Theoretically, the covariance matrix can be extracted from a test-retest dataset for the radioligand at hand. Several procedures to estimate the covariance matrix were considered as part of the simulation experiments, and the effect of the test-retest sample size was also assessed. The results are conclusive in that MLE can be used to estimate Δ and VND from PET data, avoiding the limitations associated with linear regression. The performance of LEO was, naturally, dependent on the procedure used to estimate the covariance matrix, and the test-retest sample size. Given a test-retest sample size of at least 5, but preferably 10 individuals, LEO provides higher accuracy and precision than Lassen plot in the estimation of Δ and VND. We conclude that LEO is valuable in drug occupancy studies.
Introduction Neuroimaging with Positron Emission Tomography (PET) has proven to be a valuable asset for drug development. With PET, relevant characteristics of new molecular entities can be evaluated, including level of brain exposure, target engagement (occupancy), and pharmacological activity (Gunn and Rabiner, 2017; Hargreaves et al., 2015), allowing for effective screenings already at early developmental phases. Receptor occupancy, i.e., the fraction by which a drug occupies receptors in the brain, is typically estimated from two PET measurements (one taken at baseline and the other taken after administration of the drug), using a radioligand that binds specifically to the same target. Since the drug and the PET radioligand compete for the binding sites, the difference between the two scans can be used to estimate the occupancy of the drug. Further, with well-designed studies, measurement of occupancy via PET can be used to determine the dose required for the occupancy to fall within suitable therapeutic windows (i.e., the dose range required for a clinical effect, but with limited undesirable side
effects) (Farde et al., 1988). These studies typically consist of several baseline-block scan pairs, where varying doses of the drug are administered prior to the second scan. The relationship between drug plasma concentration and the corresponding receptor occupancy can then be used to determine the drug plasma concentration required to occupy 50% of the available receptors ðKi Þ, and the maximal possible occupancy of the drug (Δmax ). The procedure with which the fractional occupancy (Δ) is estimated from PET studies, however, still deserves some attention. Binding competition between the drug and the radioligand is assumed to occur only at the specific binding sites. In addition to being specifically bound to the target of interest, however, all PET radioligands experience some degree of non-specific binding, i.e., binding to other sites than the target of interest. Also, some fraction of the radioligand is always free in tissue water. The PET system alone cannot disentangle signal originating from specifically bound radioligand from the non-displaceable (¼ free þ nonspecifically bound) radioligand. Therefore, unless the non-displaceable signal can be estimated, and in turn, the amount of specifically bound
* Corresponding author. 101 Riverside Drive, New York, NY, 10032, USA. E-mail addresses: [email protected], [email protected] (M. Schain). https://doi.org/10.1016/j.neuroimage.2018.05.017 Received 17 November 2017; Received in revised form 27 March 2018; Accepted 5 May 2018 Available online 9 May 2018 1053-8119/© 2018 Elsevier Inc. All rights reserved.
M. Schain et al.
NeuroImage 178 (2018) 255–265
radioligand can be calculated, Δ is not easily accessible. For some of the available radioligands, a brain reference region with negligible concentrations of specifically bound radioligand can be designated. For such a radioligand, Δ can easily be calculated, and they are therefore not considered in this study. For radioligands for which no reference region can be designated, estimation of Δ is not straightforward. As of today, the only established procedure for this estimation is a graphical analysis approach, in which the difference between the total distribution volumes (VT) from the two scans is plotted against that of the baseline scan, and Δ and nondisplaceable distribution volume (VND) can be obtained from linear regression. This Lassen plot was initially proposed by Lassen et al. and later revisited by Cunningham et al. (Cunningham et al., 2010; Lassen et al., 1995) There are several reasons to seek alternatives to the Lassen plot. From a theoretical standpoint, the linear regression performed in such an analysis is not strictly valid. One required assumption of regression analysis is that the independent variable is observed without measurement error; when the independent variable is measured with error, this leads to the extensively studied attenuation bias (Spearman, 1904) and has been demonstrated in PET applications (Slifstein and Laruelle, 2000). Furthermore, since baseline VT estimates are a component of both the independent and the dependent variables, the errors in the independent and dependent variables are correlated, which further violates the required assumptions in regression analysis. In this study, we present an alternative concept on how to analyze PET data acquired in blocking studies using radioligands for which no reference region exists. The general idea is to find the estimates of VND and Δ that maximize the likelihood of the observed VT values, using the same assumptions as those implicitly encoded into the Lassen plot, which are: 1) Δ is the only parameter that differs between the two scans, and 2) every region of interest (ROI) that is included in the analysis has the same value for Δ and VND. With these assumptions, joint probability density functions for the two scans can be formulated, and its maximum can be determined using standard likelihood techniques. This study is divided into 4 sections. First we provide formal mathematical derivation of the proposed algorithm, with particular emphasis on how to account for the covariance structure of the data. Second, we evaluate the algorithm's performance on simulated [11C] DASB data, including a large number of possible scenarios (varying levels of true occupancy, noise levels in the VT estimates, and procedures to estimate the covariance matrix). Third, we simulate a scenario in which estimation of Ki of a drug is sought. Finally, we test the algorithm's performance on a real dataset, in which [11C]WAY-100635 measurements were conducted before and after administration of pindolol, a partial agonist for the serotonin 1A (5HT1A) receptor. In all these scenarios, the results from the new method are compared to those obtained from the Lassen plot.
of the observed data can be expressed as
univariate Gaussian distributions, with means VS;j þ VND and pffiffiffiffiffiffi ð1 ΔÞVS;j þ VND , respectively, and standard deviation Σj;j . Further, baseline block ; VT Þ denote the joint probability density let L ðΔ; VND ; VS VT function for Vbaseline and Vblock . The analytical expression for the logT T likelihood function ℓ, neglecting the constant term that does not depend on observed data or unknown parameters, is as follows:
T
ℓ Δ; VND ; VS jVbaseline ; Vblock VND 1 VS Σ1 Vbaseline ¼ Vbaseline T T T T
VND 1 VS T VND 1 ð1 ΔÞVS Σ1 Vblock Vblock T T VND 1 ð1 ΔÞVS (3) As a function of the free parameters, equation (3) can be thought of as a surface in a kþ2 dimensional space. The coordinates at which ℓ reaches a maximum thus represent the optimal estimate for VND, Δ, and the k estimates of VS (one per ROI). Direct maximization over high dimensional space can be challenging, but in this case, it is possible to apply profile likelihood techniques and greatly reduce the dimensionality of the problem. To maximize (3) we can solve ∂V∂ S ℓ ¼ 0. After some rearrangements, we find that VS ¼
Vbaseline VND þ ð1 ΔÞ Vblock VND T T 1 þ ð1 ΔÞ2
:
(4)
For clarity purposes, we define 8 baseline > ¼ Vbaseline VND T < VS > :
Vblock ¼ S
Vblock VND T 1Δ
;
(5)
and equation (4) can be expressed as VS ¼
Vbaseline þ ð1 ΔÞ2 Vblock S S 1 þ ð1 ΔÞ2
:
(6)
Thus, given VND and Δ, there exist closed-form expressions for the maximum likelihood estimates of the other parameters, so that the loglikelihood function defined in equation (3) becomes a function of only two parameters (VND and Δ). The maximum of this function can be determined using a standard optimization routine. In our implementation, we used fminsearch. m in MATLAB R2016a. As expected, equation (6) shows that each entry of VS is replaced with the weighted average of VS estimated from the two scans, with weights determined by the occupancy. As proof of concept, Fig. 1 shows the surface spanned by equation (3) with the substitutions shown in equation b ¼ I. In addition to a global (4) in the absence of noise, and with Σ
Theory Likelihood estimation of occupancy (LEO) Let Vbaseline and Vblock denote the k 1 vectors of estimated VT values T T before and after administration of the drug in k ROIs for a given subject. According to a kinetic model with an arbitrary number of compartments, ¼ VS þ VND þ εbaseline Vbaseline T ; Vblock ¼ ð1 ΔÞVS þ VND þ εblock T
(2)
where 1 denotes a k 1 vector of ones, and Σ represents the covariance matrix of the data, whose estimation is discussed in a separate section (for now, we will treat Σ as if it is known). This means that, for ROI j, baseline block the marginal probability density functions for VT;j and VT;j are
Material and methods
Vbaseline NðVND 1 þ VS ; ΣÞ T ; Vblock NðVND 1 þ ð1 ΔÞVS ; ΣÞ T
minimum at the correct coordinates ðΔ ¼ 0:5; VND ¼ 2Þ; a curved ridge with a local minimum is observed in the parameter space, possibly reflecting that other combinations of VND and Δ may provide reasonable fits in some instances.
(1)
where Δ denotes the fractional occupancy of the drug, 0 Δ 1, VS is the distribution volume of specifically bound radioligand, and εbaseline and εblock are error terms associated with each measurement. In this model, equivalently to Lassen plot, VS varies across the k ROIs, whereas Δ and VND are constants. Assuming normality of the errors, the distribution
Estimation of covariance matrix The covariance matrix Σ in equations (2) and (3) represents the within-subject covariance matrix of the regional VT values across ROIs. This can be interpreted as the covariances across regional VT values for 256
M. Schain et al.
NeuroImage 178 (2018) 255–265
the sampled eigenvalues, and derive a procedure to non-linearly shrink the sample covariance matrix towards the population asymptote. The shrunken matrix may therefore represent a better approximation of the real, but unknown, covariance matrix than the sample covariance matrix. The procedure introduces no strong assumptions regarding the distribution of the population eigenvalues. Note that with any of the methods to estimate Σ among those outlined above (Identity, Diagonal and Shrink), the number of ROIs included in the analysis is not limited by the number of test-retest subjects. The consideration of the covariance between measurement errors in different ROIs marks a distinct difference to an earlier study that pursued a similar idea (Radua et al., 2011). In that study, the maximum likelihood formulation considered independent Gaussian noise with a fixed variance in all ROIs. In the current context, such procedure corresponds to using the identity matrix multiplied by a free parameter (σ ), whose value is estimated by the maximum likelihood estimator. However, this constant term does not affect the coordinates of the minimum of the log-likelihood function, and therefore, the previously proposed procedb ¼ I. The possibility to estimate Σ ure is in effect the same as setting Σ
Fig. 1. LEO proof-of-concept. The figure shows the surface spanned by ℓ over a range of values for Δ and VND. The figure shows ℓ rather than ℓ for visualization purposes. A minimum is observed at the correct coordinates ðΔ ¼ 0:5; VND ¼ 2Þ.
directly from test-retest data, without requiring strong assumptions of measurement noise characteristics, marks an additional distinction from the approach outlined in a recent study, which proposed to use a noise propagation model to estimate the elements of Σ (Naganawa et al., 2017).
one subject who has undergone repeated PET measurements, so that the differences in VT across scans are in large attributable to measurement errors, rather than biological variability. Since Σ contains covariances between all ROIs, it is a k k matrix of k unknown parameters (k þ parameters in the most general case). It 2 is therefore not generally feasible with two PET scans to estimate the elements of Σ in the same analysis in which we are estimating VND and Δ. Fortunately, for many radioligands, data from a test-retest experiment, in which each of several subjects are imaged twice, are available. In a testretest experiment biological variability between the test and retest scan is normally assumed to be negligible. It can be shown, that if Vtest and T Vretest denote the k N matrices, where each of the N columns contains k T regional VT values for one of the N subject included in the test-retest experiment, then the sought covariance matrix becomes Σ covðVtest T Vretest Þ=2 if N is sufficiently large (see the supplementary material for T derivation). However, unless the number of subjects included in the test-retest study is large, the covariances between ROIs may not be accurately estimated. In addition, if the number of test-retest subjects is smaller than the number of ROIs included in the analysis (i.e., N < k), the estimate of Σ will not be invertible. Most test-retest studies include 5–10 research subjects, and in such a case, we cannot count on the raw estimate of Σ being sufficiently accurate. In this study, we evaluate 3 different procedures to approximate Σ from a small-sampled test-retest study, each b with increasing complexity (estimates of Σ are denoted Σ).
Simulating data with known covariance structure It is likely that the estimation of the covariance matrix will impact the performance of LEO. To provide a “gold standard”, unattainable in practice, we will consider the situation in simulations in which the “true” covariance structure of the data is known. Let Σ denote the true, known covariance matrix for a given data set. The Cholesky factorization of Σ produces an upper triangular matrix Γ, so that Σ ¼ ΓT Γ. Let B denote a k M matrix, where each column holds k independent random numbers, sampled from Nð0; 1Þ. Then, covðΓT BÞ ¼ ΓT covðBÞΓ ¼ ΓT IΓ ¼ Σ: Let Vtrue denote a k 1 vector representing true T (i.e., noise free) VT values in k ROIs. Then, by adding the columns in ΓT B true are created, each with unique noise, and to Vtrue T , M realizations of VT with covariance matrix equal to Σ. Depending on how Vtrue is defined, T this procedure allows for simulation of both blocking and test-retest data.
Simulation of [11C]DASB blocking data In this simulation, we seek to investigate the overall performance of the various options for LEO in comparison to the Lassen plot, in terms of accuracy and precision in estimates of Δ and VND. We evaluate all 3 b (identity, diagonal, and shrink), various levels of different versions of Σ noise, a range of values for the true occupancy, and different number of b is subjects included in the hypothetical test-retest study from which the Σ obtained. A flowchart describing the simulation setup is reported in Fig. 2. This simulation is based on parameters estimated from a previously published test-retest dataset with the serotonin transporter radioligand [11C]DASB in 10 subjects (Ogden et al., 2007) (test and retest scans were was defined as the conducted on the same day). From this dataset, Vtrue T average VT across individuals and scans in k ¼ 8 ROIs that, based on previous publications, encompass a range of VT-values (Parsey et al., 2006) (anterior cingulate cortex, amygdala, dorsal putamen, insula, occipital cortex, orbitofrontal cortex, thalamus, and ventral striatum). Regional VT values were calculated using the 1-tissue compartment model, as previously suggested (Ogden et al., 2007; Ginovart et al., 2001). Further, an underlying “true” covariance matrix, Σ, was defined Vretest Þ=2, where Vtest and Vretest denote 8 10 via Σ ¼ covðVtest T T T T matrices holding regional VT values from the test and retest scans, respectively. We take this estimate of Σ to be the “true” Σ for the purposes of our simulation.
b ¼ I. Setting Σ b to equal the identity matrix corresponds to 1) Identity, Σ disregarding all differences in variance among regions and any covariance between regions. b ¼ diagðcovðVtest Vretest Þ=2Þ. Setting Σ to equal the di2) Diagonal, Σ T
T
agonal version of the estimated covariance matrix corresponds to considering distinct variances for each ROI, but setting to zero the covariances between ROIs. b obtained from non-linear shrinkage of the sample covari3) Shrink, Σ ance matrix towards the population asymptotes. Option 3 above was obtained from the work of Ledoit and Wolf, who derived a shrinkage formula for covariance matrices estimated with small sample sizes (Ledoit and Wolf, 2015). In brief, Ledoit and Wolf present a method that consistently approximates the population eigenvalues from
257
M. Schain et al.
NeuroImage 178 (2018) 255–265
above procedure. In each simulation set, 1000 instances were used (i.e., in (7), i ¼ 1; …; 1000). Parameters that varied among the different simulations were: occupancy (Δtrue ¼ 0.25, 0.50, and 0.75), noise levels (λ ¼ 1, 1.1, and 1.25, corresponding to 100%, 110% and 125% of realistic b (ranging levels), and number of test-retest subjects used to calculate Σ between 5 and 100). For each noise instance, LEO and Lassen plot were b and V d used to estimate Δ and VND (estimates denoted Δ ND ). LEO was employed with three different estimates of Σ described in 2.1.2. As a proof of concept, LEO was also employed with the Oracle covariance matrix (i.e., the matrix obtained from the measured test-retest data, which represents the true Σ in this simulation). For each simulation, and b the squared distance d between the estimated and for each choice of Σ, the dΔ ¼
true value for occupancy and VND was P1000 true P 2 b 2 ; dVND ¼ 1000 ðV true V d ΔÞ ND Þ . ND i¼1 ðΔ i¼1
calculated,
Simulation of a blocking experiment to determine Ki In this section, we simulate a scenario in which the inhibition constant Ki for a drug is sought. In this scenario, it is common that following the baseline scan, research subjects are administered different doses of the drug for the blocking scan. The baseline and blocking scans are used to estimate the Δ at each dose. With Δ estimated at different dose levels, Ki and the maximal occupancy, Δmax , can be fitted in
Fig. 2. Flow-chart describing the setup for the [11C]DASB simulation.
With Σ and Vtrue defined, the procedure outlined in the previous T section (Simulating data with known covariance structure) was used to simulate new test-retest subjects (i.e., the procedure was performed twice with the same values for Vtrue T ). These simulated test-retest subb following the procedures described in jects were then used to obtain Σ
Δ ¼ Δmax
administered before the second scan (Parsey et al., 2006). Vtrue;baseline T and Vtrue;block was defined as the average VT value in each ROI across T true and subjects. Subsequently, a Lassen plot was performed to define VND . As before, Σ, defined from the measured test-retest dataset Vtrue S described above, was used to create M noise instances with known covariance. Separate noise instances were generated for the baseline and true , Vtrue blocking data sets, ðΓT BÞbaseline and ðΓT BÞblock . Finally, with VND S , and Σ defined, the i'th realization of k simulated VT values at baseline and after block was calculated as baseline
true Vsim;baseline ¼ Vtrue þ VND 1 þ ðΓT BÞi S T;i
Vsim;block T;i
true
¼ ð1 Δ
Þ
Vtrue S
þ
true VND 1
λ
block
þ ðΓT BÞi
λ
;
(8)
where Cp is the concentration of the drug in blood plasma at the time of each of the block scans. Here, we define a hypothetical drug that binds specifically to the 5HT1A receptor, has a Ki of 25 ng/mL, and potential to occupy all receptors (Δmax ¼ 1). In our design, a total of 10 blocking experiments are conducted at varying doses of the drug. At the time of blocking scans, Cp was set to 5, 10, 20, 30, 40, 50, 60, 70, 80, and 90 ng/mL. Insertion of these numbers into equation (8) provides the “true” values for Δ that corresponds to the defined Cp. The procedure to simulate the PET data was similar to that outlined in the previous section. First, VT-values from a previous same-day testretest study (Parsey et al., 2000) in 5 individuals with the 5HT1A radioligand [11C]WAY-100635 were used to define Σ, and subsequently, to simulate a test-retest dataset. The true VT values were defined as the average VT across subjects and scans, in 4 arbitrarily selected ROIs associated with a range of values (anterior and posterior cingulate cortex, hippocampus, and insula). Note that, in the current analysis, the included number of ROIs had to be < 5, since only 5 test-retest subjects were available (with more ROIs, the “true” Σ would not have been invertible, which was required for noise generation). The true VND was set to 0.3, which is in the range of VND reported for [11C] WAY-10063515, and the noise scale factor λ was set to 1, corresponding to 100% of realistic noise. The procedure for this simulation, shown in Fig. 3, is as follows. First, 10 test-retest-subjects were simulated using Σ. Second, 10 baselineblocking scan pairs were simulated using equation (7), each with unique noise (also from Σ), and with Δ defined from equation (8) using b was calculated from the simulated testthe hypothetical Cp. Third, Σ
section Estimation of covariance matrix. Theoretically, the estimate of Σ will improve with increasing number of test-retest subjects. Therefore, the experiment was repeated with varying number of test-retest subjects, ranging from 5 to 30 in steps of 5, and then from 30 to 100 in steps of 10. Most brain PET test-retest studies include between 5 and 10 subjects, and therefore these settings provide information on the performance of the algorithm under realistic conditions. With larger testretest sample sizes, the errors introduced by a poorly estimated covariance matrix should hypothetically decrease, and the results gradually approach those obtained with a perfectly estimated (i.e., a true) covariance matrix. For this purpose, we also evaluated LEO's performance when the true covariance matrix was used, which is referred to as the “Oracle” method. This scenario represents an upper limit of LEO's performance, since it is insensitive to potential errors in the estimate of Σ. As such, the Oracle method enables us to distinguish between the effects of a poorly estimated covariance matrix and the uncertainty intrinsic to LEO. For the blocking data, we based our simulation on 8 previously published blocking experiments with [11C]DASB, where sertraline was
(
Cp Cp þ Ki
retest measurements, using only Shrink (the other procedures to estib were omitted in results since previous simulation was conclusive mate Σ b Fourth, VND and in that Shrink was the preferred procedure to obtain Σ). Δ were estimated both using LEO and Lassen plot. Finally, using the known Cp and the estimated values for Δ, Ki and Δmax were estimated using equation (8). Repeating this procedure 1000 times results in 1000 estimates of Ki and Δmax obtained with Lassen plot, and, 1000 estimates for LEO. The distributions of Ki and Δmax were visualized with histograms, and the agreement were quantified by calculating the squared difference
(7)
where ðΓT BÞi denote the i'th column in the noise matrix, Δtrue represents the true occupancy, and λ is a scale factor for the noise. A total of 27 sets of simulations were performed, each following the 258
M. Schain et al.
d dKi ¼
between the estimated and the true P1000 true P max 2 bi Þ2 ; dΔmax ¼ 1000 ðΔmax;true Δd ðK K Þ . i i¼1 i¼1
NeuroImage 178 (2018) 255–265
Δ was also calculated using the Lassen plot and LEO, and the results were compared to those estimated using the reference region (average across ROIs). In contrast to the simulations above, the true b was esticovariance matrix Σ is not known. When applying LEO, Σ
values,
mated from the [11C]WAY-100635 test-retest dataset in 5 individuals that was used in section Simulation of a blocking experiment to determine b was only estimated with the Shrink method, since preKi14. Here, Σ
Application to real data In this final section, we apply Lassen plot and LEO to a previously published dataset (Martinez et al., 2001), consisting of [11C] WAY-100635 measurements at baseline, and after administration of pindolol, a partial 5HT1A receptor agonist. In this dataset, arterial input functions were acquired, allowing for estimation of regional VT. In addition, for [11C]WAY-100635, the cerebellar white matter has been shown to be a suitable reference region, allowing for estimation of VND. To estimate VND, a 1-tissue compartment model was fitted to the time-activity curve for the cerebellar white matter. To estimate regional VT values, a 2-tissue compartment model was applied, with K1/k2 constrained to the estimate of VND, as previously suggested (Parsey et al., 2000). Estimates of VT and VND enables estimation of regional binding potential (Parsey et al., 2005). Here, we used BPP (¼ VT-VND) as outcome measure (Parsey et al., 2000). Using the reference region, the occupancy of pindolol was estimated via Δ ¼ ðBPbaseline P Þ=BPbaseline in 7 ROIs (amygdala, posterior cingulate cortex, BPblock P P hippocampus, occipital cortex, orbitofrontal cortex, parietal cortex, and parahippocampal gyrus). These ROIs were selected since they, after preliminary analysis of BPP, displayed similar occupancy within each subject (<10% difference from across-ROI mean).
vious analyses showed that this method is preferable over the simpler procedures to estimate Σ. Occupancies obtained using the reference region were compared to those obtained from Lassen plot and LEO, via regression analysis.
Implementation notes In the presence of noise, both LEO and the Lassen plot can “fail”, i.e., provide highly unrealistic estimates. To avoid this, constraints on the free parameters were enforced in LEO. By definition, 0 Δ 1, so the occupancy was constrained to be in this interval by performing an uncon Δ strained minimization over log 1Δ rather than over Δ. With a similar reasoning, a non-negativity constraint was enforced for VND, by minimizing over logðVND Þ. Although the Lassen plot has an implicit upper bound on the occupancy, it is generally executed without explicit constraints on the free parameters. Thus, the mere presence of parameter constraints in LEO has the potential to result in lower scores of dΔ and dVND , without actual improved accuracy. To avoid this, all occupancy and VND values obtained from the Lassen plot that were <0 were then set to 0 before dΔ and dVND were calculated. Also, the “failure rate” was defined as the number of simulation instances for which a method provided an estimate of VND that was larger than the maximal VT for that simulated subject. All simulations that failed according to this definition were excluded before calculation of dΔ and dVND . Results Simulation of [11C]DASB blocking data Results from simulation 1 are shown in Figs. 4–7. Figs. 4 and 5 report the sum of squared distances d between true and estimated values for Δ and VND for various test-retest sample sizes (thus, a low value reflects a good performance). LEO using the true underlying covariance matrix (i.e., the “Oracle” method) always resulted in the lowest squared distances. Although this is interesting from a theoretical point of view, the “Oracle” method is not possible in practice, as it requires knowledge of the true underlying covariance structure of the data. Of interest, with increasing test-retest sample sizes, using non-linear shrinkage to estimate the covariance matrix provided results that asymptotically approached those obtained with the Oracle method. Considering commonly used test-retest sample sizes (i.e., 10 or 5 subjects), the results indicate that LEO performed with covariance matrices estimated with either the Shrink method, or diagonal matrices, is preferable over Lassen plot, as these methods consistently resulted in lower scores of dΔ and dVND . The difference between LEO and Lassen plot b was less evident when the identity matrix was used in place of Σ. To allow further visualization of some of the results shown in Figs. 4 and 5, the distributions of the estimates of Δ and VND obtained with 10 test-retest subjects are reported in Figs. 6 and 7, respectively. In these violin plots, the relative width and color within each violin represent the density of the respective distributions. The black dashed line corresponds to the true value in each plot. The peaks around 0 seen for Lassen plot and LEO executed with the identity and diagonal matrix for 25% occupancy in Fig. 6 results from the lower bound constraint; all instances for which the methods would provide a negative occupancy
Fig. 3. Procedure for the simulation of a blocking experiment to determine Ki.
259
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 4. Sum of squared error in estimates of occupancy (dΔ ) for various noise levels, occupancy levels, and test-retest sample sizes used to estimate the covariance matrix Σ. The lowest error was consistently observed with the Shrink method, which asymptotically approached results obtained with the true Σ (i.e., Oracle). Shrink also resulted in smaller error than Lassen plot even for small sample sizes.
are floored at 0 (see Implementation notes). Although no clear bias was observed in occupancy estimated with the Lassen plot, there was a considerable spread. This spread was reduced with LEO when executed with Shrink approach to estimate Σ. In contrast to what was observed for occupancy, all methods resulted in a systematic overestimation of VND (Fig. 7). The largest bias was observed at 25% occupancy, with Lassen plot and LEO executed with the identity matrix in place of the covariance matrix. For higher occupancy levels, the bias was reduced. In all cases, LEO executed with the Shrink approach to approximate Σ resulted in lower bias and lower spread than the other procedures available in reality. In supplementary figure S1, the true and approximated covariance matrices are visualized for varying number of test-retest subjects. Inspection of these matrices reveals that Σ is not approximately diagonal, but rather, many off-diagonal elements are not negligible. Further, figure S2 in the supplementary material shows associations of the diagonal elements of Σ with average and standard deviation of the VT values (calculated from the baseline data), and with ROI volumes. These associations were r2 ¼ 0.98 for average VT values, r2 ¼ 0.93 for their standard deviations, and r2 ¼ 0.60 for ROI volumes (Figure S2, supplementary material).
Lastly, Table S1 in the supplementary material reports the rate by which each of the evaluated methods “fails” according to the criterion described in Implementation notes. For all methods, highest failure rates were observed at low occupancy levels (for instance, at 100% noise, Lassen plot failed in 5% of the instances at 25% occupancy, and never at 75% occupancy). Shrink, regardless of test-retest sample size, consistently resulted in the lowest failure rate (with 10 test-retest subjects and 100% noise, the failure rate was 1.4% at 25% occupancy). The highest b failure rate was observed when the identity matrix was used in place of Σ, whereas Lassen plot and LEO using a Diagonal matrix resulted in similar failure rates.
Simulation of [11C]WAY-100635 blocking data to estimate Ki From the previous experiment, it was evident that using the Shrink approach to estimate the covariance matrix in LEO provided better results than when the diagonal or identity matrix was used. For this purpose, these simplified versions were omitted from the remaining experiments. The “Oracle” method was also excluded since it is not available in real data applications.
260
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 5. Sum of squared error in estimates of VND (/1000) (dVND ) for various noise levels, occupancy levels, and test-retest sample sizes used to estimate the covariance matrix Σ. The lowest error was consistently observed in correspondence of the Shrink method, which asymptotically approached results obtained with the true Σ (i.e., Oracle). Shrink also resulted in smaller error than Lassen plot even for small sample sizes.
were r2 ¼ 0.99, slope ¼ 1.02, and intercept ¼ 0.04. The Lassen plots carried out to estimate the occupancies shown in Fig. 9 are presented in Figure S3 in the supplementary material.
Following the procedure outlined in Fig. 3 resulted in 1000 datasets, each with 10 estimated values for Δ corresponding to the defined drug plasma levels. The estimated Δ were then used in equation (8) to estimate Ki and Δmax . With Ki and Δmax estimated, the occupancy plot describing the relationship between drug plasma levels and receptor occupancy were obtained. The top panel in Fig. 8 display these relationships with occupancy estimated with LEO and the Lassen plot, respectively, and the bottom panels show the distributions of Ki and Δmax . LEO resulted in a lower spread in the estimates of Ki and Δmax compared with the Lassen plot. The squared distances for the two parameters were dKi ¼ 211:5; and dΔmax ¼ 0:034 for Lassen plot, and dKi ¼ 59:0 and dΔmax ¼ 0:010 for LEO.
Discussion Graphical procedures for the analysis of PET data are attractive, due to their simplicity to implement and interpret (Logan et al., 1990; Patlak et al., 1983). The Lassen plot is one of the most commonly used approaches to estimate occupancy using radioligands for which no suitable reference region exist, in spite of it being not strictly valid due to correlated errors and errors in the predictor variable. To avoid the assumptions implicit in linear regression, we here consider maximum-likelihood (ML). ML is a well-suited method for the problem at hand, as it relies on a solid theoretical framework assuring that the estimates of the free parameters of interest are approximately unbiased. In addition, the ML estimator is consistent, meaning that it is capable of determining the true parameter value with arbitrarily high precision given a sufficient number of measurements. In the derivation of a ML estimator for occupancy and VND, it is important not to introduce additional assumptions regarding the underlying theory of the pharmacokinetic properties of the radioligand or the drug. For this purpose, we formulated the likelihood function based
Application to real data As a final experiment, we applied LEO and Lassen plot to a previously published dataset, in which 9 subjects had [11C]WAY-100635 PET measurements obtained before and after administration of pindolol. The occupancy calculated using LEO (with Shrink) or the Lassen Plot was compared to that obtained from using BPP. The scatter plot comparing these occupancy levels is shown in Fig. 9. The linear regression analysis resulted for Lassen plot in r2 ¼ 0.83, slope ¼ 1.33, and intercept ¼ 0.20. For LEO, the corresponding values 261
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 6. Violin plots showing the distribution of occupancy estimated with Lassen plot and LEO, assuming 10 test-retest subjects. The relative width and color in each violin represents the density, and the dashed black line corresponds to the true value. No evident bias was observed for any method, but LEO executed with Shrink resulted in the lowest spread. The additional yellow fields seen at 25% occupancy result from non-negativity constraints (see section Implementation notes).
Although a perfectly acquired test-retest dataset may provide a theoretical possibility to calculate the covariance matrix, there are practical limitations. The accuracy of the estimated covariance matrix depends on the sample size of the test-retest dataset. From a test-retest sample size of 5–10 individuals, simply calculating Σ using the definition for covariance is not recommendable, as the covariances will be poorly estimated. In addition, if the number of ROIs included in the analysis is larger than the number of test-retest subjects, Σ will not be invertible, which disables application of LEO. Rather, alternative methods to estimate Σ need to be undertaken. Note that, with such methods to estimate Σ (for instance, using the diagonal matrix or nonlinear shrinkage), the number of ROIs included in the analysis is not constrained by the test-retest sample size. Here, we evaluated three such methods, and the results show that estimation of the covariance matrix is of importance. First, if equal variances and zero covariances were assumed (i.e., the identity matrix was used in place of Σ), LEO resulted in similar spread and bias as that observed for the Lassen plot. Second, if different within ROI variances were considered but covariances between ROIs were still neglected (i.e., using the diagonal matrix), a somewhat improved estimation of Δ and VND was observed. The performance improved even further when the full covariance matrix was estimated via non-linear shrinkage. Finally, in the hypothetical situation that the true covariance matrix was completely known, VND and occupancy were estimated with high accuracy and precision. This latter option can be thought of as the upper bound of the LEO's performance, which is very unlikely to ever be reached in a real application. For real clinical investigations, the results from this study indicate that the non-linear shrink method to estimate the covariance matrix based on 10 test-retest subjects is favorable in comparison the Lassen plot.
on the standard assumptions in compartmental modeling: each molecule of the radioligand that has penetrated the brain is, at any given time point, either specifically bound to the target, or non-displaceable (Innis et al., 2007). Further, the error term was assumed to follow a normal distribution. This latter assumption is in large supported by the fact that LEO takes as input the estimated VT-values, which are typically derived from non-linear regression. According to likelihood theory, it then follows that the errors in VT are (at least asymptotically) normally distributed. However, the error associated with an estimate of VT may come from multiple sources, for instance statistical noise in the emission and transmission data, error in the image reconstruction, and measurement error of the arterial input function. All of these error sources have the possibility to influence different brain regions similarly. It is therefore not reasonable to assume that the errors in two separate brain regions from the same subject are independent. To allow for dependency in errors across brain regions, the covariance structure across brain regions was considered in the derivation of LEO. In a blocking experiment, each subjects VT values are compared at baseline and after administration of the drug, and the fractional occupancy, which is treated as a free parameter, is the only parameter that is assumed to differ between the scans (i.e., VS and VND are assumed equal between baseline and block). The covariance matrix should therefore only contain covariance between measurement errors, without the influence of biological variability that is present across subjects. We here propose to approximate this information using test-retest datasets for the radioligand, since the difference between a test and a retest scan in one individual should in large reflect only measurement errors. However, under certain circumstances, some biological variability between test and retest can be expected due to, for instance, diurnal variation (Matheson et al., 2015). 262
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 7. Violin plots showing the distribution of VND estimated with Lassen plot and LEO, assuming 10 test-retest subjects. The relative width and color in each violin represents the density, and the dashed black line corresponds to the true value. For 25% occupancy, an evident overestimation in VND was observed for all methods (except the Oracle method), and the bias was reduced at higher occupancy levels. Of the methods that are typically available in real studies, LEO executed with Shrink provided the smallest bias and spread.
at hand. That LEO explicitly uses test-retest information possibly enables estimation of Δ at lower values. Last, we investigated the agreement between estimates of occupancy obtained using a reference region to those obtained using Lassen plot and LEO. As seen in Fig. 9, occupancy estimated with LEO was on group level in closer agreement with that estimated using a reference region than estimates obtained with Lassen plot. These results should however be treated with some caution, since the sample is rather small, and it is possible that, at least in part, the improved performance seen with LEO is driven by the parameter constraints. When comparing LEO and Lassen plot using real data, no particular attention was given to estimates of VND. For [11C]WAY-100635, VND is very low, typically around 0.3 when estimated using the VT in the cerebellar white matter (Parsey et al., 2005). For 5 out of the 9 subjects, Lassen plot resulted in negative values of VND (figure S3 in supplementary material), whereas with LEO these values were positive (possibly as a result from the non-negativity constraint), but very close to 0. From these observations, we considered [11C]WAY-100635 VND estimated with either method to be biologically not meaningful. Regardless of whether Lassen plot or LEO is used, it is important to be aware that both algorithms assume each ROI included in the analysis to have the same value for VND and occupancy. It is not obvious how to verify this assumption using observed data. We note that a violation of this assumption can result in an unpredictable bias that may confound findings. In the Lassen plot, a non-linear relationship or clear outliers may give hints as to whether certain ROIs should be excluded from the analysis. We suggest that when estimating parameters using LEO, that also visually examining the Lassen plot can
An interesting question is how well the covariance matrix can be estimated in a scenario where test-retest data are not available. Supplementary figure S2 shows that, using only single scan data, the diagonal elements of Σ can be very well predicted using VT values and ROI volumes, indicating that single scan data can likely be used to estimate within-ROI variances. However, as seen in supplementary figure S1, Σ is not approximately diagonal – many off-diagonal elements have nonnegligible values. Further research is needed to understand how to better approximate off-diagonal elements of Σ either using single scan data, or, as recently proposed, modeling the noise characteristics of the PET data (Naganawa et al., 2017). To assess the impact of using LEO in real life applications, we simulated a scenario in which the properties of an unknown new drug is sought, with the underlying hypothesis that the estimation Ki and Δmax will improve as a result of improved accuracy and precision of Δ. From the results in Simulation 1, the largest difference between Lassen plot and LEO was observed at low or intermediate occupancy values, where LEO provided higher accuracy and precision in the parameter estimates. Poorly estimated occupancies at low drug plasma levels (typically corresponding to low occupancies) impacts the entire fit of equation (8) and may thus produce uncertainties in the estimates of Ki, regardless of the quality of the estimates at higher plasma levels. We therefore hypothesized that, with LEO, the estimation of Ki will improve as a result of improved accuracy and precision of Δ, in particular at low drug plasma levels. The results presented in Fig. 8 confirmed this hypothesis, where a lower spread across noise instances is observed in the occupancy curves (top panel). In general, it is likely that the lower limit of a detectable occupancy is determined by the test-retest variability of the radioligand
263
M. Schain et al.
NeuroImage 178 (2018) 255–265
Fig. 8. Results from the simulated [11C]WAY-100635 experiment. Top panel show the relationship between receptor occupancy and plasma concentration obtained with LEO and Lassen plot (each line represents one simulation). Bottom panels show the distributions of estimates of Ki and Δmax . In each plot, the black line corresponds to the true value.
Funding This study was supported by the Paul Janssen Fellowship in Translational Neuroscience. Acknowledgements We thank Prof. Michael Wolf and Dr. Olivier Ledoit at the Department of Economics, University of Zurich, Switzerland, for helpful discussions regarding estimation of covariance matrices from small datasets. Appendix A. Supplementary data Supplementary data related to this article can be found at https://doi. org/10.1016/j.neuroimage.2018.05.017. References Cunningham, V.J., Rabiner, E.A., Slifstein, M., Laruelle, M., Gunn, R.N., 2010. Measuring drug occupancy in the absence of a reference region: the Lassen plot re-visited. J. Cereb. Blood Flow. Metab. 30, 46–50. Farde, L., Wiesel, F.A., Halldin, C., Sedvall, G., 1988. Central D2-dopamine receptor occupancy in schizophrenic patients treated with antipsychotic drugs. Arch. Gen. Psychiatry 45, 71–76. Ginovart, N., Wilson, A.A., Meyer, J.H., Hussey, D., Houle, S., 2001. Positron emission tomography quantification of [(11)C]-DASB binding to the human serotonin transporter: modeling strategies. J. Cereb. Blood Flow. Metab. 21, 1342–1353. Gunn, R.N., Rabiner, E.A., 2017. Imaging in central nervous system drug discovery. Semin. Nucl. Med. 47, 89–98. Hargreaves, R.J., Hoppin, J., Sevigny, J., et al., 2015. Optimizing central nervous system drug development using molecular imaging. Clin. Pharmacol. Ther. 98, 47–60. 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. 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.
Fig. 9. Scatter plot between occupancies obtained using a reference region to calculate BPp, and those obtained from VT values using Lassen plot and LEO. Σ was estimated using the Shrink method applied to VT values from 5 testretest subjects.
provide a good complementary analysis and potentially reveal problems with the data. Finally, to facilitate the application of LEO for other researchers; we freely distribute MATLAB code for doing the calculations. The code and instructions are available at https://github.com/martinschain/LEO. 264
M. Schain et al.
NeuroImage 178 (2018) 255–265 humans: comparison of arterial and reference tisssue input functions. J. Cereb. Blood Flow. Metab. 20, 1111–1133. Parsey, R.V., Arango, V., Olvet, D.M., Oquendo, M.A., Van Heertum, R.L., John Mann, J., 2005. Regional heterogeneity of 5-HT1A receptors in human cerebellum as assessed by positron emission tomography. J. Cereb. Blood Flow. Metab. 25, 785–793. Parsey, R.V., Kent, J.M., Oquendo, M.A., et al., 2006. Acute occupancy of brain serotonin transporter by sertraline as measured by [11C]DASB and positron emission tomography. Biol. Psychiatry 59, 821–828. Patlak, C.S., Blasberg, R.G., Fenstermacher, J.D., 1983. Graphical evaluation of blood-tobrain transfer constants from multiple-time uptake data. J. Cereb. Blood Flow. Metab. 3, 1–7. Radua, J., Bullich, S., Lopez, N., Catafau, A.M., 2011. Restricted maximum likelihood estimation of PET neuroreceptor occupancy in the absence of a reference region. Med. Phys. 38, 2558–2562. Slifstein, M., Laruelle, M., 2000. Effects of statistical noise on graphic analysis of PET neuroreceptor studies. J. Nucl. Med. 41, 2083–2088. Spearman, C., 1904. The proof and measurement of association between two things. Am. J. Psychol. 15, 72–101.
Ledoit, O., Wolf, M., 2015. Spectrum estimation: a unified framework for covariance matrix estimation and PCA in large dimensions. J. Multivar. Anal. 139, 360–384. Logan, J., Fowler, J.S., Volkow, N.D., et al., 1990. Graphical analysis of reversible radioligand binding from time-activity measurements applied to [N-11C-methyl](-)-cocaine PET studies in human subjects. J. Cereb. Blood Flow. Metab. 10, 740–747. Martinez, D., Hwang, D., Mawlawi, O., et al., 2001. Differential occupancy of somatodendritic and postsynaptic 5HT(1A) receptors by pindolol: a dose-occupancy study with [11C]WAY 100635 and positron emission tomography in humans. Neuropsychopharmacology 24, 209–229. Matheson, G.J., Schain, M., Almeida, R., et al., 2015. Diurnal and seasonal variation of the brain serotonin system in healthy male subjects. Neuroimage 112, 225–231. Naganawa, M., Gallezot, J.D., Rossano, S., Carson, R.E., 2017. Quantitative PET imaging in drug development: estimation of target occupancy. Bull. Math. Biol. Ogden, R.T., Ojha, A., Erlandsson, K., Oquendo, M.A., Mann, J.J., Parsey, R.V., 2007. In vivo quantification of serotonin transporters using [(11)C]DASB and positron emission tomography in humans: modeling considerations. J. Cereb. Blood Flow. Metab. 27, 205–217. Parsey, R.V., Slifstein, M., Hwang, D.R., et al., 2000. Validation and reproducibility of measurement of 5-HT1A receptor parameters with [carbonyl-11C]WAY-100635 in
265