Target and Tissue Selectivity Prediction by Integrated Mechanistic Pharmacokinetic-Target Binding and Quantitative Structure Activity Modeling

Selectivity is an important attribute of effective and safe drugs, and prediction of in vivo target and tissue selectivity would likely improve drug development success rates. However, a lack of understanding of the underlying (pharmacological) mechanisms and availability of directly applicable predictive methods complicates the prediction of selectivity. We explore the value of combining physiologically based pharmacokinetic (PBPK) modeling with quantitative structure-activity relationship (QSAR) modeling to predict the influence of the target dissociation constant (K D) and the target dissociation rate constant on target and tissue selectivity. The K D values of CB1 ligands in the ChEMBL database are predicted by QSAR random forest (RF) modeling for the CB1 receptor and known off-targets (TRPV1, mGlu5, 5-HT1a). Of these CB1 ligands, rimonabant, CP-55940, and Δ8-tetrahydrocanabinol, one of the active ingredients of cannabis, were selected for simulations of target occupancy for CB1, TRPV1, mGlu5, and 5-HT1a in three brain regions, to illustrate the principles of the combined PBPK-QSAR modeling. Our combined PBPK and target binding modeling demonstrated that the optimal values of the K D and k off for target and tissue selectivity were dependent on target concentration and tissue distribution kinetics. Interestingly, if the target concentration is high and the perfusion of the target site is low, the optimal K D value is often not the lowest K D value, suggesting that optimization towards high drug-target affinity can decrease the benefit-risk ratio. The presented integrative structure-pharmacokinetic-pharmacodynamic modeling provides an improved understanding of tissue and target selectivity.


INTRODUCTION
In the development of new therapeutics, the balance between safety and efficacy is critical for success. The selectivity of a new drug is therefore an important attribute, as selective compounds are less likely to mediate side effects (1). On the other hand, targeting multiple targets simultaneously is increasingly considered as a valuable option to exert sufficient effect on a complex biological system (2,3). Regardless of the desired degree of selectivity, understanding and prediction of the target binding to multiple targets in multiple tissues is essential for the optimization of pharmacotherapy. In this article, we differentiate between two types of selectivity: target selectivity and tissue selectivity. Target selectivity is defined as a difference in target binding to different receptors, and tissue selectivity is defined as a difference in target binding to the same target in different tissues. Additionally, a distinction is made between equilibrium selectivity and kinetic selectivity. Equilibrium selectivity refers to differential target binding while target binding is in equilibrium with the free drug concentration around the target. This equilibrium binding is described for single step target binding without target turnover according to Eq. 1, in which K D is the dissociation constant, [L] is the unbound drug concentration, [R] is the unbound target concentration, [LR] is the bound drug-target complex concentration, k off is the first order target dissociation rate constant, and k on is the second order target association rate constant: Equilibrium target selectivity is thus driven by differential K D values for the different targets. Kinetic selectivity, however, is often used to describe that different k off values for different targets can lead to a difference in the duration of target occupancy (4). Unfortunately, kinetic selectivity is not clearly defined and is also used to describe a difference in in vitro k off values (5). However, differential k off values do not always result in a differential duration of target occupancy in vivo since the plasma and local pharmacokinetics can also be rate-limiting for the duration of target occupancy (6,7). To disambiguate, we define Bin vitro kinetic selectivity^here as a difference in in vitro k off values and we define Bin vivo kinetic selectivity^as a difference in target occupancy between different targets that cannot be explained by a difference in equilibrium binding (e.g., a difference in K D or target site distribution).
A previous study that analyzed a minimal mechanistic model for drug elimination, tissue distribution, and target binding showed that an increase in drug-target affinity decreases the chance of observing in vivo kinetic selectivity, especially for slow tissue distribution and a high target concentration (7). On that basis, it is expected that the optimal K D for target and tissue selectivity is dependent on the target concentration, tissue distribution kinetics, and binding kinetics. This contrasts with the current practice of drug discovery and development which often aims at a minimal value for the K D and k off and a maximal ratio to the off-target K D and k off value if selectivity is concerned. The minimal mechanistic model that was analyzed in the study of de Witte et al. (7) did not consider (i) the effects of slow distribution of a drug into tissues where no target binding takes place nor (ii) the limiting role that blood flow can have on tissue distribution. In order to capture the influence of these pharmacokinetic mechanisms, physiologically based pharmacokinetic (PBPK) models can be used. In these models, a distinction is made between system-specific properties and drug-specific properties. In this type of analysis, the values of system-specific parameters such as blood flows and volumes for each organ are based on the physiological literature data, while the values of drug-specific parameters, such as partition coefficients and protein binding are often based on in vitro data or on quantitative structure activity relationships (QSARs) (8). As such, these models allow the prediction of plasma and tissue unbound drug concentrations. The influence of drug-target binding on free drug concentrations has been described frequently with target-mediated drug disposition (TMDD) models (9). The combination of PBPK and TMDD modeling has been reported in the literature previously but is not generally used in selectivity optimization (10)(11)(12)(13). To obtain the drug-specific properties that determine the values of the partitioning parameters in PBPK models, either experimental data for each individual drug or quantitative structure-activity relationships (QSAR) are required. These QSARs enable the prediction of partitioning parameters from the molecular structure. While these QSARs are often used in PBPK modeling to predict non-specific tissue distribution parameters, the prediction of specific target binding parameters is currently not incorporated in PBPK modeling, based on the assumption that the amount of drug bound to its biological target is negligible relative to the total amount of drug in the body (14)(15)(16)(17).
QSAR models may be either regression or classification models which predict a response variable from a set of predictor values. In regression models, these predictor values are related to a continuous response variable (e.g., a K D value), while in classification models the predictor values relate to a categorical variable (e.g., labeled Bactive^or Binactive^). The predictor values represent the molecular structure and molecular properties, and the response variable is an activity value, such as the K D in the case of affinity.
Machine learning methods such as support vector machines (SVMs), decision trees such as random forests (RFs) and deep neural networks (DNNs) are generally used to obtain a predictive learning model (18)(19)(20). The training of these models is based on prior data, which means that their performance is greatly dependent on data quality and availability. A suitable database for bioactivity data is available in the ChEMBL, which can be used to obtain predictive QSAR models (21,22).
Integration of drug-target binding prediction and pharmacokinetic modeling allows for the prediction of the selectivity profile for a given ligand directly from its molecular structure. As such, this modeling approach may provide information on a ligand's efficacy and safety in vivo during the early stages of drug development. This is especially relevant in systems that contain off-targets or targets that are also expressed in organs were no drug effect is desired. An example of the latter system is the cannabinoid system, of which the cannabinoid receptor CB1 is a major component. The CB1 receptor is widely expressed throughout the body but mainly found in the brain where it mediates a broad range of effects in health and disease (23,24). Many off-targets have been identified for CB1 ligands, including the vanilloid receptor TRPV1, the metabotropic glutamate receptor mGlu5, and the serotonin receptor 5-HT1a (25,26). Activity at these receptors, predominantly in the brain, may amplify or counteract effects at the CB1 receptor. TRPV1, for example, has been suggested to have an effect opposite of that of CB1 in anxiety and depression, which are common side effects observed for CB1 antagonists, and mGlu5 is a major player in the GABA-system, which is the target system for CB1mediated therapies in Parkinson's disease (27)(28)(29). The mechanisms underlying functional in vivo selectivity are diverse and complex, but computational elucidation of offtarget affinities and their integration in combined PBPK-TMDD modeling could help to identify safety concerns early in drug discovery and development, and to design the optimal properties of new drugs.
This article describes an approach towards the development of an integrative predictive modeling for drug selectivity. Firstly, the main determinants of in vivo equilibrium and kinetic selectivity are identified by minimal PBPK-TMDD modeling and simulation. Secondly, the development and validation of a Random Forest-based QSAR (QSAR-RF) model for the prediction of K D values is described. Lastly, an example of the use of predicted K D values in PBPK-TMDD modeling is provided for the combined in vivo target and tissue selectivity of rimonabant, a prototype antagonist at the CB1 receptor.

Software
All simulations were performed in RStudio Version 1.0.136 coupled to R version 3.4.0 (30,31). Physicochemical property prediction and QSAR modeling were performed in Pipeline Pilot version 2016 (32).

Pharmacological Models
Three PBPK-TMDD models were developed: a minimal PBPK-TMDD model for simulation of target selectivity (model I,

Parameters
Models I and II. All physiological values of the systemspecific parameters were obtained from literature (33)(34)(35)(36)(37)(38). The heart was used as a reference organ for the determination of distribution into and out of the binding tissue. An overview of all model parameters is supplied in Supplemental 1. Blood flows were summed for all tissues that are lumped and converted to rate constants using the tissue volumes.
Model III. All physiological values of the system-specific parameters were obtained from literature (33)(34)(35)(36)(37)(38). Target site distribution in the brain was not characterized by the blood flow but by the average effective flow through the target site as obtained from literature values of clearance out of the brain extra-cellular fluid due to the fluid flow as estimated for nine drugs (39). The conversion of these values as well as an overview of all parameters are supplied in Supplemental 1.
Receptor densities of CB1, mGlu5, TRPV1, and 5-HT1a in the cerebellum, hypothalamus, and frontal cortex were obtained from the literature for all four receptors, except the receptor concentration of mGlu5 in the hypothalamus and 5-HT1a in the cerebellum, which were not reported in the literature (40)(41)(42)(43)(44). The mGlu5 receptor concentration in the hypothalamus was filled in with the averages of the other brain regions since differences between brain regions for the other receptors did not differ drastically. The 5-HT1a receptor concentration in the cerebellum was set to the low value of 0.01 nM as it was reported to be unidentifiable (43). Receptor concentrations in rats and humans were used interchangeably since no complete set of receptor densities could be obtained for either rats or humans. Values found in literature have shown to differ no more than tenfold (41,45). TRPV1 concentrations were given in nanogram/milligram lysate and converted to picomole/milligram protein by linear conversion. For this, the receptor concentration in nanogram/ milligram lysate and femtomole/milligram protein in the hypothalamus as reported in the literature was used (42,46). The receptor density in the hypothalamus in femtomole/ milligram was divided by the receptor density in nanogram/ milligram lysate, and the resulting coefficient was used to transform the receptor density in nanogram/milligram lysate of the cerebellum, hypothalamus, and frontal cortex to the corresponding receptor density in femtomole/milligram. CB1 and TRPV1 concentrations in picomole/milligram were then converted to nanomolar using a conservative (i.e., the lowest published) estimate of protein concentration in wet tissue of 100 mg/mL from literature (47)(48)(49). An overview of the target concentrations is presented in Table I. An overview of the conversions and all target concentrations can be found in Supplemental 2.
Tissue-blood partition coefficients were calculated according to Poulin & Theil 2000 (Eq. (2)) (50), and water and lipid fractions were obtained from Poulin & Krishnan1995 (51). The required physicochemical parameters (logP, logSo) (52) were determined in Pipeline Pilot. An overview of all parameters is supplied in Supplemental 3. where: Pt:b predicted value of the tissue-blood partition coefficient So the solubility of the ligand in n-octanol (mol*m-3) Sw the solubility of the ligand in water (mol*m-3) Nb the neutral lipid content of blood (as fraction of blood volume) Nt the neutral lipid content of the tissue (as fraction of tissue volume) Pb the phospholipid content of blood (as fraction of blood volume) Pt the phospholipid content of the tissue (as fraction of tissue volume) Wb the water content of blood (as fraction of blood volume) Wt the water content of the tissue (as fraction of tissue volume)

Simulations
Model I. Model I was used to investigate the influence of K D , target concentration (R tot ), and k off on in vivo target selectivity. To this end, four different simulations (a, b, c, d) were performed. In all four simulations, the k off at the first target (R1) was set to 0.01 h −1 and the k off at the second target (R2) was set to 10 h −1 while both the K D and R tot were the same for both targets. An overview of the parameter values that were varied in these simulations can be found in Table II. An overview of all other parameters can be found in Supplemental 1.
Model II. This model was used to perform simulations to investigate the influence of K D , target concentration (R tot ), and tissue distribution (k in ) on in vivo tissue selectivity. To this end, four different simulations were performed for a k in value of 8.6 h −1 (fast tissue distribution) and for a k in value of 0.86 h −1 (slow tissue distribution). An overview of the variable parameter values can be found in Table II. An overview of all other parameters can be found in Supplemental 1.
Model III. Simulations were performed for rimonabant, Δ 8 tetrahydrocannabinol (Δ 8 THC), and CP-55940 in a minimal PBPK-TMDD model (Fig. III). The K D at the selected targets (CB1, mGlu5, TRPV1, and 5-HT1a) was predicted by a QSAR per target model trained on the complete pChEMBL dataset per target. A fast dissociation from the receptor was assumed for all compounds by setting the k off value to 10 h −1 at all receptors. Simulations were performed for a time span of 7 days during which a dose was administered every 24 h.
In order to investigate the influence of increasing drugtarget affinity without a change in equilibrium selectivity, additional simulations were performed in which the ratio between the different K D values for the different receptors was kept the same while adjusting the absolute K D values by a factor of 10 and 100. Simulations were performed for a time span of 7 days with dosing once every 24 h. The dose was scaled for the K D to obtain similar equilibrium occupancies in all simulations.

QSAR
A Random Forest QSAR per target model was developed using the Random Forest package from CRAN (53).

Data Selection
Bioactivity data from ChEMBL22 where used for the development of the QSAR model (54). High-quality data was selected by setting assay confidence at 9 and requiring an assigned pChEMBL value for all data points (22). This means that a direct single protein target is assigned to the ligand. PubChem database data and potential duplicates were excluded from the dataset. Bioactivity data from ChEMBL was limited to four different constants: K D , K i , IC 50 , and EC 50 . It has been shown previously that K i and IC 50 can be combined for modeling (55). In order to check if these constants could be used interchangeably, a statistical analysis of their pChEMBL values was performed. In this analysis, the mean, standard deviation (SD), and median and median absolute deviation (MAD) were analyzed within and between all four constants. An overview of all results is provided in Supplemental 4. Since from this analysis it could be concluded that the deviation between pChEMBL values between K D , nbt non-binding tissue, c central compartment, bt binding tissue, el eliminating tissue and K i do not differ significantly from the deviation within the K D dataset, both K D and K i values were used in the model development.
The molecular structure of the ligands was extracted from the molfile and physicochemical properties, and FCFP_6 circular fingerprints were calculated in Pipeline Pilot (56). The FCFP_6 fingerprints were then converted to 768 feature properties for use in model training. Selection was performed based on the relative frequency of substructures per target, where the optimal frequency was close to being present in 50% of the ligands.
The complete dataset was split into a training set (70%) and validation set (30%). This split was performed seven times, each time with a different seed (111, 222, …, 777) in order to create seven different datasets. In this way, the model training and validation could be performed seven times, allowing for reproducibility analysis of the model performance results.

Training
For each target, a Random Forest model consisting of 500 trees was trained using the seven different training sets. The models were trained on a predefined set of properties consisting of log(P), molecular weight, number of proton donors, number of proton acceptors, number of rotatable bonds, number of atoms, number of rings, number of aromatic rings, molecular solubility, molecular surface area, molecular polar surface area, and the 768 FCFP_6 fingerprint properties that describe the molecular structure in more detail.

Validation
The model performance was validated internally and externally using the corresponding validation dataset per seed, as described above. Internal validation was performed by an out-of-bag (OOB) estimate and presented as the average R 2 regression coefficient and the rootmean-squared error (RMSE) (57). The OOB estimate method uses subsamples from the training dataset to determine the mean prediction error of the RF model. The RMSE is a value that measures the average magnitude of the error and is presented by the same unit as the dependent variable, which in this case is the pChEMBL value (-log K D /K i in M). External validation

Model I
The simulations in Fig. IV show in vivo kinetic target selectivity in all simulations, due to a difference in the k off value for target 1 and target 2. However, the extent of the observed selectivity is dependent on the K D value and target concentration. Given that optimization is often performed towards lower k off values, the target at which k off is 0.01 h −1 is considered as the desired therapeutic target. Initially, target selectivity for the off-target is observed, but this selectivity reverses to selectivity for the therapeutic target over time in all simulations, except in Fig. IVb, where the K D is low and the target concentration is high. As it would be unlikely in drug development to develop two drugs with a 1000-fold different binding kinetics but the same K D value, we also performed these simulations with 100-fold different binding kinetics and 10-fold different K D values as presented in Supplemental 5. In these figures, the same trend is observed, although the residual selectivity is higher due to the difference in K D values. In addition, the same simulations as for Fig. IV were performed with repeated dosing until steady-state was reached. In these simulations (Fig. S3), a similar reduction in selectivity is observed as in Fig. IV, only to a smaller extent.
In summary, we observed that the combination of a high target concentration and a low K D value decreases the in vivo kinetic target selectivity.

Model II
For the simulations presented in Fig. V, no difference in k off or K D values between the two target sites could be included, since the ligand binds to the same target and the differences in target occupancy arise due to a difference in the target concentration. No selectivity is observed for the higher K D values (10 and 1 nM), and only marginal selectivity is observed for lower K D values (0.1 and 0.01 nM). Assuming that the target concentration in the target tissue is higher than the target concentration in the   off-target tissue, the lowest simulated K D values showed selectivity in the first 12 h to the off-target tissue after which selectivity for the target tissue is observed (Fig. Vd).
Marginal selectivity for the target tissue is observed for a K D value of 0.1 nM (Fig. Vc). Taken together, this means that the K D and receptor concentrations influence the extent of in vivo tissue selectivity. In addition, the same simulations as for Fig. V were performed with repeated dosing until steady-state was reached. In these simulations (Fig. S4), the reduction in peak occupancy for tissue 1 compared to tissue 2, as observed in Fig. V, is reduced and the average occupancy becomes approximately equal in steady-state. The simulations in Fig. V were performed for fast tissue distribution based on the reported blood flow of wellperfused organs in the human body (37). Figure VI shows the simulation results for slower tissue distribution, representing limited perfusion of the target site (e.g., in a synaptic cleft) or the presence of diffusion barriers (e.g., for intracellular or CNS targets). In these figures, the same patterns are observed as for fast tissue distribution, but the observed selectivity is greater and the affinity for maximal selectivity for the target-rich tissue is lower. Similarly as for

QSAR-RF
From the simulations described above, it follows that there is an optimal K D for both tissue selectivity and target selectivity. To facilitate the optimization of the K D , we aimed to predict the K D value from the molecular structure with predictive QSAR modeling. In this study, a QSAR-RF model was developed. The results of the internal and external validation are given in These values indicate good model performance, since the error in public data is around 0.44 for pK i data. Moreover, based on this error, it has been shown that the theoretical maximal achievable R 2 value then becomes 0.81 for the perfect model (58)(59)(60). A full overview of the results is supplied in Supplemental 4.

Model III
To reflect a drug discovery/candidate selection scenario, the developed QSAR model was used to predict the affinity of the molecules rimonabant, Δ 8 tetrahydrocannabinol (Δ 8 THC), and CP-55940 for the four selected receptors (CB1, TRPV1, mGlu5, 5-HT1a, Fig. VIII). These K D values were then used to predict the selectivity over different brain regions (cerebellum, hypothalamus, and the frontal cortex). The results of these simulations are given in Fig. VIII. For the target occupancy of Δ 8 THC, the compound with the lowest CB1 affinity, no selectivity is observed between brain regions. The target occupancy for the higher affinity compounds rimonabant and CP-55490 shows a slower increase of target occupancy in the brain regions with the highest target concentrations, the cerebellum and frontal cortex compared to the hypothalamus. The difference in target occupancy between the brain regions is similar for all targets, which results in a change in target selectivity across brain regions. Two days after the start of rimonabant dosing, for example, the simulated target occupancy at TRPV1 in the hypothalamus is similar to the CB1 target occupancy in the cerebellum and frontal cortex.
For CP-55940, it takes more than 7 days to reach the maximal occupancies in the cerebellum and frontal cortex, while this delay would be even more extensive for lower doses. It should be noted that equilibrium selectivity (i.e., the difference in K D values for the different receptors) is different for the compounds in Fig. VIII. To obtain a better view of the role of the value of the K D as such, rather than the K D ratio between targets, the simulations for rimonabant were repeated with the same K D ratio between targets and tenfold increased and decreased K D values, as shown in Supplemental 6, Fig. S6. Additionally, to explore the influence of error propagation from the QSAR model into model III, simulations were performed for the lowest and highest K D value within the RMSE-based K D (±) prediction range as shown in Supplemental 6, Fig. S7. Summarizing the results, it is consistently found that the selectivity profile changes drastically over time while this would not be expected based on plasma concentrations and K D values alone.

DISCUSSION
In this study, we have shown that the integration of target binding and PBPK modeling demonstrates the importance of target concentrations, target site distribution kinetics, and the K D and k off for both in vivo target selectivity and tissue selectivity. We observe that a low K D , in combination with a high target concentration, decreased the kinetic target selectivity, indicating that the duration of target occupancy is less sensitive to the value of k off for these conditions. Moreover, we find that an increasing K D can both increase and decrease tissue selectivity, dependent on the target concentration and tissue distribution. This influence of the K D on tissue selectivity is especially relevant in the first days after a change in dosing regimen as the difference in average target occupancy between a relatively high and low K D compound disappears upon reaching steady-state. The demonstrated mechanistic modeling can thus be instrumental to find an optimal K D value for a specific target/therapeutic area. To utlize this approach most effectively, our QSAR model to predict CB1 and off-target K D values can be used to predict tissue and target selectivity directly from the molecular structure. Using this combination of models, our simulations for the CB1 ligands confirm that lower K D values for all targets can decrease the CB1 and brain region selectivity significantly during the first days of treatment.
Our results suggest that optimization towards high drugtarget affinity and slow drug-target dissociation, as is commonly performed within the current drug development paradigm, may not result in the most selective compounds. While this study demonstrates the influence of target concentrations on the target occupancy in different tissues, the influence of target concentrations on the occupancyresponse relationship has previously been described as the driving factor for tissue selectivity of partial agonists (61)(62)(63). Also, the influence of target expression and target accesibility on target saturation and its duration has been described in quantitative terms before (64). Moreover, the relevance of the k off for transient tissue selectivity has been demonstrated before in a PBPK model for inhaled drugs (65). For the development of more selective drugs, target concentrations of both the intended target and off-targets as well as distribution to the target tissue/target site should be taken into consideration. In this respect, it is important to consider that distribution to the target site is not only dependent on distribution into the target tissue, but also on the localization of the target within this tissue (e.g., in the bloodstream or intracellularly). Moreover, factors such as target concentrations and tissue distribution may be altered in a disease state, Predicted K D values of Δ 8 -THC, rimonabant, and CP-55940 at the CB1, 5-HT1a, mGlu5, and TRPV1 receptor were used in these simulations. k off values were assumed to be 10 h −1 . The drug dose was scaled linearly for the CB1 K D value and was administered every 24 h. R tot,cer,CB1 = 527 nM, R tot,cer,mGlu5 = 5.1 nM, R tot,cer,TRPV1 = 19 nM, R tot,cer,5-HT1a = 0.01, R tot,hyp,CB1 = 248 nM, R tot,hyp,mGlu5 = 16 nM, R tot,hyp,TRPV1 = 13 nM, R tot,hyp,5-HT1a = 2.37, R tot,fc,CB1 = 529 nM, R tot,fc,mGlu5 = 25 nM, R tot,fc,TRPV1 = 12 nM, R tot,fc,5-HT1a = 1.7 nM which is important for the translation from healthy volunteers to patients (24,(66)(67)(68). Finally, it should be considered that there is an increased interest towards allosteric modulation in CNS drug discovery due to the potential benefits with regard to selectivity and side effects (69). However, it has also been shown that allosteric modulators display different physicochemical and efficacy (K i versus ligand efficacy) profiles compared to orthosteric ligands (70). These parameters can be included in the modeling approach for future studies.
The methods described in this study provide valuable insights for drugs in later stages of the drug development process. Especially for tissue selectivity, but also for target selectivity, the intended duration of therapy needs to be taken into account to define optimal drug properties and appropriate study design. Drug therapies that are meant for shortterm dosing or to have a quick initiation of the drug effect, such as pain killers for acute pain, anti-inflammatory and antiinfectuous drugs, and anesthetics, are especially influenced by the impact of target concentrations and affinity on tissue selectivity. For drug therapies that are meant for chronic dosing, the influence of the target concentration and the drugtarget affinity should be taken into account for evaluation of drug candidates in the initial phase after the first and the last drug dose. The selectivity profiles in Fig. VIII, for example, would result in underestimation of CB1 selectivity in (pre)clinical studies, if only the first 7 days were studied. This might lead to the unnecessary discontinuation of the development of valuable drug candidates. Moreover, the slowly increasing target occupancy for high-affinity drugs such as CP-55940 might lead to a clinically undesired delay and unfavorable selectivity between the initiation of treatment and the onset of the therapeutic effect. While for CP-55940 and the current dosing regime the delay is in the order of a week which is probably acceptable in the clinic, this delay will extend to less acceptable periods in the order of months for drugs with even higher affinities (e.g., biologics) or for drugs that require lower doses. This can potentially be mitigated by a higher dose (i.e., a loading dose), which can be lowered as soon as steady-state occupancy is reached. Since monitoring of occupancy levels in the clinic is hardly feasible, this would required in-depth knowledge of the mechanisms and predicted occupancy profile as described in this study. Moreover, it should be noted that the target occupancy will decline only slowly after discontinuation of treatment and that it might take several days or even weeks for these high-affinity drugs before the target occupancy is back to insignificant levels. This could be counteracted in the clinic by administration of a competitive antagonist or agonist to displace the drug from the receptor and enhance the clearance out of the target binding tissue. In Fig. VIII, we simulated the target occupancy in different brain areas, where the CB1 receptor is expressed in high concentrations. However, the CB1 receptor is expressed at much lower concentrations in other organs, which are also involved in the influence of CB1 antagonism on metabolic disorders (71,72). Our results suggest that more selectivity for peripheral CB1 receptors and less psychological side effects could be achieved by dosing for a few days and subsequently stopping the drug dosing until the central CB1 occupancy is lowered to minimal levels, after which the sequence can be repeated. The clinical benefit of such a regimen is entirely hypothetical, but could be worth exploring.
The simulations in this study were all based on physiological parameter values as obtained from PBPK models and target concentration literature. However, additional assumptions were sometimes necessary. For the simulations in Fig. VI, the tissue distribution of the drug was not based on the blood flow through well-perfused organs, as for the other figures, but we assumed a delayed distribution due to, for example, limited diffusion into a synaptic cleft or the cytosol. The magnitude of this delay is compound and target specific, and this assumption will thus only hold for a limited number of compounds. Secondly, the simulation in Fig. VIII assumed fast binding kinetics as the actual binding kinetics of rimonabant have been reported to be complex and therefore hard to accurately determine in in vitro studies (73,74). The assumption of fast binding kinetics is supported by the short dissociation half-life as reported by Packeu et al. (75). Additionally, this assumption will be valid for any drug for which the binding kinetics are not rate-limiting compared to the pharmacokinetics, but slower binding kinetics could change the outcome of the simulations, as shown in Fig. IV and in previous studies (6,7). Thirdly, a number of assumptions concerning (interspecies) translatability of target densities were made in order to obtain useful target densities for the simulations in Fig. VIII. In general, the quality of absolute tissue-specific target concentration data, rather than relative expression values, might be limited. This is illustrated by the large deviations between experimental tissue density results found in the literature between PET-studies and tissue Bno wash^assay experiments (41,76). In addition, basic input parameters such as brain region volumes already come with uncertainty and variability (77). Furthermore, the limited amount of information on target site distribution for the simulations in Fig. VIII limits the predictive value of these simulations. Moreover, protein binding is not incorporated in the simulations underlying Fig. VIII and the target binding in our model is driven by unbound + nonspecifically bound tissue concentrations, while for many drugs the unbound tissue concentrations drive the drug-target binding. However, it should be noted that the dose in Fig. VIII is adjusted to get a desired occupancy. Inclusion of protein binding and the free concentration as the driving factor for target binding can be achieved by multiplying the free + nonspecifically bound drug concentration with the fraction unbound. To get the same target occupancy as presented in our simulations, the dose would need to be adjusted accordingly. However, in most cases, the kinetics and selectivity of target binding will not be affected by protein binding and tissue partitioning, especially not compared to our simulations for Fig. VIII since tissue partitioning is assumed to be the same in all brain areas in these simulations. These simulations should therefore be considered as a prediction of the relevant parameters for combined target and tissue selectivity for a realistic set of target concentrations and K D values, rather than a precise prediction of target occupancy values for the simulated CB1 ligands. One of the most striking findings in our study is that increasing the K D in drug development can both increase and decrease the target and tissue selectivity. This demonstrates the relevance of target concentrations and tissue distribution, and the valuable role of mechanistic modeling.
The prediction error that is observed for the K D predictions of the developed QSAR model introduces an extra level of uncertainty into the overall reliability of the selectivity predictions. The largest RMSE value in this study was found for the mGlu5 QSAR, with an average value of 0.8. This value relates to the deviation of the predictions from the actual pChEMBL value, and has the same unit as the dependent variable, which in this case is the -log K D . This uncertainty is therefore carried on into the pharmacological simulations. From the simulations performed with the highest and lowest value within the K D prediction range of rimonabant, it can be concluded that this propagation of error does influence the observed selectivity profile. This error is limited to the extent of selectivity and the distribution across brain regions during the first 1 to 4 days. However, part of this error is already present in the public data that was used to train our QSAR model, in which a larger standard deviation is found compared to the rimonabant predictions at the CB1 receptor from the QSAR model (Supplemental 4, Fig. SI). Additionally, having the ability to predict the selectivity profile in the earliest stages of drug disovery justifies the use of predictions with significant uncertainty. Moreover, both the overrepresentation and underrepresentation of structural features or scaffolds in the ChEMBL database might decrease the predictive power for new compounds that do not share these structural features.
Although the predictive value of the presented models is limited by the assumptions we made, the presented insight into the influence of the target concentration and tissue distribution kinetics is in line with the previous analysis of more simple models with only one target and one tissue (7). Moreover, the relevance of incorporating target binding in PBPK models for the accurate prediction of tissue concentrations has been demonstrated before (13). The basic principle behind the role of the K D and target concentration on the duration of occupancy is the high concentration of the target compared to unbound drug concentrations at the target site. This leads to a depletion of the free drug concentration during drug-target association after initiation of drug administration and a replenishment of the free drug concentration during drug-target dissociation after termination of the drug administration. The depletion of free drug concentrations can be local and can be hard to observe from systemic plasma concentrations, as demonstrated for the simulations of CP-55940 in Fig. S8. This is mainly relevant for drugs with a high K D and target concentration and at a target occupancy that is not completely saturated. If this target occupancy is increasing, drug-target association will deplete the unbound target site concentration, and if the occupancy is decreasing, drugtarget dissociation will increase the unbound target site concentration compared to plasma concentrations.
In summary, the information presented in this study provides new insights into the mechanisms underlying in vivo target and tissue selectivity, specifically in relation to drugtarget affinity, target concentration, tissue and target site distribution, as well as binding kinetics. The study provides situations in which selectivity is expected to occur, which may aid as a lead towards creating ligands with the desired selectivity profile. Additionally, the demonstrated integration of mechanistic modeling and machine learning could enable the incorporation of these insights in the earliest phases of drug discovery. The need for this model-based selectivity optimization is especially valuable for therapeutic areas in which an optimal target or tissue selectivity profile is difficult to obtain (e.g., in oncology) and might be less valuable for therapeutic areas were selectivity is less challenging and the traditional minimization of the K D is desired (e.g., for antibiotic/antiviral targets that are not expressed in human cells).

CONCLUSIONS
Simulations performed in semi-physiological pharmacological models with target binding revealed an important role for the target concentration and tissue distribution, next to the K D and k off values, in determining the extent of selectivity. Interestingly, it was observed that the optimal selectivity is not observed for the drug that displays the highest drug-target affinity when assuming that the desired target concentrations are high and the desired binding kinetics are slow. Additionally, it was observed that kinetic selectivity is unlikely when the target concentrations and the drug-target affinity are high, while tissue selectivity is first increased and then decreased for increasing target concentrations and drug-target affinities. The context-dependent optimum of drug-target affinity in determining the extent of selectivity demonstrates the value of K D prediction for drug development. Taken together, this study demonstrates the potential of integrative predictive modeling in providing improved strategies to optimize drug candidates for maximal in vivo selectivity.