A digital human brain phantom (image size = 121 × 145 × 121, 1 mm isometric voxels) was used to create occupancy images with regional variation in occupancy, corresponding to different plasma concentrations. The brain phantom data was released under the Creative Commons Attribution-NonCommercial license (CC BY-NC) with no end date. Original MRI scans are from OASIS (https://www.oasis-brains.org/). Labelings were provided by Neuromorphometrics, Inc. (http://Neuromorphometrics.com/) under academic subscription. Concentration–response curves were generated according to an Emax model to generate 10 ideal occupancy images for the 10 different plasma concentrations of drug.

$$\text=}_}* \frac}_}$$

(1)

In Eq. 1, Occ is the occupancy in every voxel of the occupancy image; EC50 is the predefined EC50 value at every voxel of EC50 image; Occmax is the maximum occupancy at every voxel; and C is the plasma concentration. The 10 different values of C were selected between 3 – 279 ng/mL based on an existing human data set, and Occmax was set to 0.85 [16].

To generate regional variation in occupancy image, the ground truth EC50 image (shown in Fig. 1) was created with an EC50 = 25 ng/mL in caudate and EC50 = 50 ng/mL in putamen. The rest of the brain was assigned an EC50 value between 6 and 10 ng/mL consistent with a whole brain average value from human data shown in a previous study [16].

Ground truth EC50 image showing “hot spots” of EC50 in putamen and caudate regions in coronal, axial, and sagittal views

We have shown previously that noise in occupancy images is a function of occupancy as follows [14].

$$_}=-0.22*\text+0.25$$

(2)

In Eq. 2, \(_}\) is the standard deviation in occupancy, and Occ is the occupancy [14]. This noise was applied based on a normal random distribution at every voxel to all idealized occupancy images to the generate noisy occupancy images.

To add correlation between the voxels to our noisy occupancy images, we applied a Gaussian filter with a kernel size of 3 × 3 × 3. Multiple phantoms with different amounts of voxel correlation were generated using different Gaussian standard deviations (\(_}\)) for the Gaussian filter. The occupancy images (shown in Fig. 2) generated with \(_}\) = 0.5 voxels best represented the smoothness observed in occupancy images produced from real human data [14].

Simulated noisy smoothed occupancy images at different plasma concentrations. The occupancy noise model (Eq. 2) was applied to ideal EC50 images (Fig. 1) and smoothed by a Gaussian filter. White text for each image shows the plasma concentration used to generate the occupancy image from the true EC50 phantom

SLIC algorithm was used to combine multiple voxels of the occupancy images into super-voxels (clusters) [17]. SLIC is an adaptation of k-means for super-pixel generation. It uses a smaller search area in its distance calculation which is faster than other k-means algorithms. It also uses a weighted average of spatial and feature-space distances which can be used to emphasize one of the distances over the other.

The SLIC algorithm was first introduced by Achanta et al. [17] for its faster speed, greater memory efficiency and better adherence to boundaries compared to other k-means clustering algorithms. Later, a modified version of SLIC, ‘SLICR’, was introduced to incorporate temporal features from 2D dynamic computed tomography myocardial perfusion imaging [18]. In SLIC-Occ, we modified SLIC by introducing a distance in feature-space, feature-space refers to occupancy and spatial distances are calculated in 3D. The total distance measure to be minimized was calculated as:

$$D= \sqrt_}^+_}}\right)}^* ^}$$

(3)

In the Eq. 3, m controls the weighting of spatial distance over the feature distance.

$$_}= \sqrt_- _\right)}^+ _- _\right)}^+ _- _\right)}^}$$

(4)

where (xc, yc, zc) is the coordinate of the center of the cluster c and (xi, yi, zi) is the coordinate of the voxel i which is to be assigned to a cluster. Feature-space is made up of \(V\) occupancy images corresponding to different plasma concentrations. Distance in feature-space is:

$$_}=\sqrt_^_^-_^\right)}^}$$

(5)

where, \(_^\) is the occupancy value at the center of cluster c corresponding to plasma concentration value v and, \(_^\) is the occupancy value of the voxel i corresponding to plasma concentration value v. The search area for every voxel was defined as 2S x 2S x 2S, where S is defined as:

In Eq. 6, N is the total number voxels in each 3D occupancy image and the image is divided into K clusters.

There are two parameters (i.e., m and K) that need to be optimized for SLIC-Occ clustering algorithm. The number of initial clusters, K, will determine the approximate size of a cluster, N/K, in terms of voxels. While the shape of each cluster will be determined by the value of m; the larger the m the more regular the clusters.

We created multiple simulations with different combinations of m and K to investigate their effects on accuracy and precision of EC50 images in our clustering algorithm. The choice of m and K was made to reduce the CV(EC50) while maintaining accuracy of EC50.

SLIC-Occ was used to segment 10 occupancy images, corresponding to 10 different plasma drug concentrations, into super-voxels (K clusters). Average occupancy of all the voxels within each cluster was used as the occupancy value for the corresponding cluster. Two versions of the Emax model (Eq. 1), were used to fit the occupancy data. In version 1, Occmax was fixed (1-parameter model); in the 2-parameter version Occmax and EC50 were estimated simultaneously.

The corrected Akaike information criterion (AICc) was calculated for both (1-parameter and 2-parameter) Emax model fits at every concentration–response curve (i.e., every cluster) as [14, 19]:

$$\text=2p+n*\text\left(\frac}\right)+ \frac^+2p}$$

(7)

where p is the number of estimated parameters in the model, n is the number of data points being fitted, and SSE is the sum of squared errors. The model with lower AICc was selected as the best model.

Parametric images were constructed with the parameter estimate of the best model for each cluster. In other words, the final parametric images generated are a combination of 1- and 2-parameter fits depending on which model was selected based on AICc for each cluster. Using the cluster-level best-fit, EC50 and CV(EC50) images were generated. The coefficient of variation for an estimated parameter was defined as:

where \(\mu\) is the parameter estimate, and \(_\) is the standard deviation of the parameter estimate calculated by:

$$_}}^=\text\left(^*J\right)}^*\left(\frac^*R}\right)\right)$$

(9)

where J is the Jacobian matrix of Emax fit at the solution, R is the vector of residuals, n is the number of fitted data points, and p is the number of estimated parameters.

The distinct image regions in the ideal EC50 image were used to select the same voxels in the estimated EC50 image for bias calculation. Accuracy was calculated as percent bias as:

$$\text= \frac_- \text }_} }_}*100$$

(10)

The two versions of the Emax model (1-parameter and 2-parameter), were used to fit the noisy occupancy data at every voxel. As with cluster-level estimates, voxel-level parametric images were constructed with the parameter estimate of the best model at each voxel, using the AICc to determine the best model fit. EC50, and CV(EC50) images for voxel-level estimation were generated to compare with the results from cluster-level images.

SLIC-Occ algorithm was applied to human occupancy data (image size = 121 × 145 × 121, 1 mm isometric voxels) that were published in previous studies [14, 16]. A detailed description of the PET acquisition has been published [16]. In short, 5 healthy subjects underwent 3 scans each at the Yale PET center for 2 h on a ECAT EXACT HR + scanner (Siemens Medical Systems, Knoxville, TN, USA) after injection with 570 ± 141 MBq (injected mass: 2.7 ± 1.3 µg) of 11C-flumazenil, a nonselective GABAA tracer. One of the scans was performed at the baseline (no drug administration), and two were acquired after oral administration of a single acute dose of the \(\alpha 1\)-, \(\alpha 2\)-, \(\alpha 3\)-selective GABA positive allosteric modulator, CVL-865 (also known as PF-06372865). The drug dose was either 10 mg (n = 3) or 65 mg (n = 2). The post-drug scans were acquired at approximately 1.5 h and 24 h after administration of the drug. The plasma concentration of CVL-865 was measured at three different time points during the scan and averaged [16]. The occupancy versus drug concentration curves were generated for each voxel in the previous study [13, 14]. Using SLIC-Occ, the occupancy images were segmented into clusters. EC50, and CV(EC50) images were generated and compared with voxel-level EC50 and CV(EC50) images.