Surely you are joking, Mr Docking!

In the wake of recent COVID-19 pandemics scientists around the world rushed to deliver numerous CADD (Computer-Aided Drug Discovery) methods and tools that could be reliably used to discover novel drug candidates against the SARS-CoV-2 virus. With that, there emerged a trend of a significant democratization of CADD that contributed to the rapid development of various COVID-19 drug candidates currently undergoing different stages of validation. On the other hand, this democratization also inadvertently led to the surge rapidly performed molecular docking studies to nominate multiple scores of novel drug candidates supported by computational arguments only. Albeit driven by best intentions, most of such studies also did not follow best practices in the field that require experience and expertise learned through multiple rigorously designed benchmarking studies and rigorous experimental validation. In this Viewpoint we reflect on recent disbalance between small number of rigorous and comprehensive studies and the proliferation of purely computational studies enabled by the ease of docking software availability. We further elaborate on the hyped oversale of CADD methods' ability to rapidly yield viable drug candidates and reiterate the critical importance of rigor and adherence to the best practices of CADD in view of recent emergence of AI and Big Data in the field.

Estimation of Maximum Recommended Therapeutic Dose Using Predicted Promiscuity and Potency.

We report a simple model that predicts the maximum recommended therapeutic dose (MRTD) of small molecule drugs based on an assessment of likely protein-drug interactions. Previously, we reported methods for computational estimation of drug promiscuity and potency. We used these concepts to build a linear model derived from 238 small molecular drugs to predict MRTD. We applied this model successfully to predict MRTDs for 16 nonsteroidal antiinflammatory drugs (NSAIDs) and 14 antiretroviral drugs. Of note, based on the estimated promiscuity of low-dose drugs (and active chemicals), we identified 83 proteins as "high-risk off-targets" (HROTs) that are often associated with low doses; the evaluation of interactions with HROTs may be useful during early phases of drug discovery. Our model helps explain the MRTD for drugs with severe adverse reactions caused by interactions with HROTs.

Drug repurposing: far beyond new targets for old drugs.

Repurposing drugs requires finding novel therapeutic indications compared to the ones for which they were already approved. This is an increasingly utilized strategy for finding novel medicines, one that capitalizes on previous investments while derisking clinical activities. This approach is of interest primarily because we continue to face significant gaps in the drug-target interactions matrix and to accumulate safety and efficacy data during clinical studies. Collecting and making publicly available as much data as possible on the target profile of drugs offer opportunities for drug repurposing, but may limit the commercial applications by patent applications. Certain clinical applications may be more feasible for repurposing than others because of marked differences in side effect tolerance. Other factors that ought to be considered when assessing drug repurposing opportunities include relevance to the disease in question and the intellectual property landscape. These activities go far beyond the identification of new targets for old drugs.

Linking pharmacology to clinical reports: cyclobenzaprine and its possible association with serotonin syndrome.

The link between cyclobenzaprine (Flexeril) administration and serotonin syndrome (SS) is subject to debate. Establishing such a connection is difficult because of the limited number of case reports available and the almost complete ignorance of its preclinical pharmacology. In this context, evidence is provided here that cyclobenzaprine blocks the serotonin and norepinephrine transporters and binds to another set of five serotonin receptors. SS should be considered when indicative signs occur in the context of cyclobenzaprine use.

Solid-phase synthesis of libraries generated from a 4-phenyl-2-carboxy-piperazine scaffold.

Strategies for finding novel structures of therapeutical interest are discussed. The rationale for the selection of the two scaffolds N4-(m-aminophenyl)-piperazine-2-carboxylic acid E and N4-(o-aminophenyl)-piperazine-2-carboxylic F is described. The synthesis of the appropriate precursors to scaffold E and F and their use in solid-phase chemistry are described. A 160-member library was produced combining these novel piperazine scaffolds with eight sulfonyl chlorides/acid chlorides and 10 amines. The compound library prepared was analyzed using LC-MS, showing the expected base peak in all wells at an average purity of 82%.

Rapid estimation of hydrophobicity for virtual combinatorial library analysis.

Novel NPH (Nonpolar, Polar, Hydrogen) descriptors for rapid estimation of hydrophobicity, amenable for filtering extremely large virtual combinatorial libraries (VCL) are proposed, based on atom counts: P_at, the sum of polar atoms (sum of oxygen, nitrogen, phosphorus and sulfur); NP_at, the sum of non-polar atoms (the sum of carbons and halogens minus P_at); SHDA (the sum of hydrogen bond donors and acceptors). In combination with molecular weight, the following related parameters are defined: MWP_at (the "polar" molecular weight); MWNP_at (the "nonpolar" molecular weight); and MWSHDA (the "hydrogen bonding" molecular weight). The NPH descriptors provide moderate-to-good predictive PLS models when external prediction is evaluated against measured logP values (q2 pred > 0.5, n = 7954) or against calculated logP values (q2 pred > 0.6, n = 18 991). Related to hydrophobicity, the NPH descriptors are intended for fast analyses of extremely large VCLs (10(6)-10(12) compounds), even before the enumeration of reactants into products occurs.

MTD-PLS: a PLS-based variant of the MTD method. A 3D-QSAR analysis of receptor affinities for a series of halogenated dibenzoxin and biphenyl derivatives.

MTD-PLS, the Partial Least Squares (PLS) variant of the Minimum Topological Difference (MTD) method is described. In MTD-PLS, molecules are characterised not only by the occupancy or nonoccupancy of the hypermolecular vertices (as in classical MTD), but also by additional descriptors for each vertex: fragmental van der Waals volumes, fragmental hydrophobicities, partial atomic charges, etc. This method was applied to a series of 73 polyhalogenated derivatives of dibenzo-p-dioxine, dibenzofuran and biphenyl (induction of aryl hydrocarbon hydrolase and affinities to rat cytosolic receptor), previously studied by MTD. The separation of steric, hydrophobic, and electrostatic effects was achieved retranslating from the latent variable space into a linear combination of the initial structural variables. The MTD-PLS method yields more detailed results compared to classical MTD, indicating the importance of electrostatic effects at some substituent positions.

Is there a difference between leads and drugs? A historical perspective.

To be considered for further development, lead structures should display the following properties: (1) simple chemical features, amenable for chemistry optimization; (2) membership to an established SAR series; (3) favorable patent situation; and (4) good absorption, distribution, metabolism, and excretion (ADME) properties. There are two distinct categories of leads: those that lack any therapeutic use (i.e., "pure" leads), and those that are marketed drugs themselves but have been altered to yield novel drugs. We have previously analyzed the design of leadlike combinatorial libraries starting from 18 lead and drug pairs of structures (S. J. Teague et al. Angew. Chem., Int. Ed. Engl. 1999, 38, 3743-3748). Here, we report results based on an extended dataset of 96 lead-drug pairs, of which 62 are lead structures that are not marketed as drugs, and 75 are drugs that are not presumably used as leads. We examined the following properties: MW (molecular weight), CMR (the calculated molecular refractivity), RNG (the number of rings), RTB (the number of rotatable bonds), the number of hydrogen bond donors (HDO) and acceptors (HAC), the calculated logarithm of the n-octanol/water partition (CLogP), the calculated logarithm of the distribution coefficient at pH 7.4 (LogD(74)), the Daylight-fingerprint druglike score (DFPS), and the property and pharmacophore features score (PPFS). The following differences were observed between the medians of drugs and leads: DeltaMW = 69; DeltaCMR = 1.8; DeltaRNG = DeltaHAC =1; DeltaRTB = 2; DeltaCLogP = 0.43; DeltaLogD(74) = 0.97; DeltaHDO = 0; DeltaDFPS = 0.15; DeltaPPFS = 0.12. Lead structures exhibit, on the average, less molecular complexity (less MW, less number of rings and rotatable bonds), are less hydrophobic (lower CLogP and LogD(74)), and less druglike (lower druglike scores). These findings indicate that the process of optimizing a lead into a drug results in more complex structures. This information should be used in the design of novel combinatorial libraries that are aimed at lead discovery.

Cheminformatics: a tool for decision-makers in drug discovery.

Cheminformatics is a tool that aims at facilitating the decision-making process across various preclinical stages of drug discovery. Access to biological and chemical data, but not the data themselves, is an integral part of cheminformatics. Emerging tools that allow storage of, and access to, chemical, structural-chemical and biological information are only now beginning to reach maturity. Recent advances in cheminformatics include virtual library analysis without enumeration and novel methods to investigate global chemical similarity and diversity voids. The most important task for cheminformatics is to constantly reevaluate itself and its utility in the area of drug discovery, in order to provide probabilistic, rather than categorical predictions.

Chemography: the art of navigating in chemical space.

Combinatorial chemistry needs focused molecular diversity applied to the druglike chemical space (drugspace). A drugspace map can be obtained by systematically applying the same conventions when examining the chemical space, in a manner similar to the Mercator convention in geography: Rules are equivalent to dimensions (e.g., longitude and latitude), while structures are equivalent to objects (e.g., cities and countries). Selected rules include size, lipophilicity, polarizability, charge, flexibility, rigidity, and hydrogen bond capacity. For these, extreme values were set, e.g., maximum molecular weight 1500, calculated negative logarithm of the octanol/water partition between -10 and 20, and up to 30 nonterminal rotatable bonds. Only S, N, O, P, and halogens were considered as elements besides C and H. Selected objects include a set of "satellite" structures and a set of representative drugs ("core" structures). Satellites, intentionally placed outside drugspace, have extreme values in one or several of the desired properties, while containing druglike chemical fragments. ChemGPS (chemical global positioning system) is a tool that combines these predefined rules and objects to provide a global drugspace map. The ChemGPS drugspace map coordinates are t-scores extracted via principal component analysis (PCA) from 72 descriptors that evaluate the above-mentioned rules on a total set of 423 satellite and core structures. Global ChemGPS scores describe well the latent structures extracted with PCA for a set of 8599 monocarboxylates, a set of 45 heteroaromatic compounds, and for 87 alpha-amino acids. ChemGPS positions novel structures in drugspace via PCA-score prediction, providing a unique mapping device for the druglike chemical space. ChemGPS scores are comparable across a large number of chemicals and do not change as new structures are predicted, making this tool a well-suited reference system for comparing multiple libraries and for keeping track of previously explored regions of the chemical space.

Property distribution of drug-related chemical databases.

The process of compound selection and prioritization is crucial for both combinatorial chemistry (CBC) and high throughput screening (HTS). Compound libraries have to be screened for unwanted chemical structures, as well as for unwanted chemical properties. Property extrema can be eliminated by using property filters, in accordance with their actual distribution. Property distribution was examined in the following compound databases: MACCS-II Drug Data Report (MDDR), Current Patents Fast-alert, Comprehensive Medicinal Chemistry, Physician Desk Reference, New Chemical Entities, and the Available Chemical Directory (ACD). The ACDF and MDDRF subsets were created by removing reactive functionalities from the ACD and MDDR databases, respectively. The ACDF subset was further filtered by keeping only molecules with a 'drug-like' score [Ajay et al., J. Med. Chem., 41 (1998) 3314; Sadowski and Kubinyi, J. Med. Chem., 41 (1998) 3325] below 0.8. The following properties were examined: molecular weight (MW), the calculated octanol/water partition coefficient (CLOGP), the number of rotatable (RTB) and rigid bonds (RGB), the number of rings (RNG), and the number of hydrogen bond donors (HDO) and acceptors (HAC). Of these, MW and CLOGP follow a Gaussian distribution, whereas all other descriptors have an asymmetric (truncated Gaussian) distribution. Four out of five compounds in ACDF and MDDRF pass the 'rule of 5' test, a probability scheme that estimates oral absorption proposed by Lipinski et al. [Adv. Drug Deliv. Rev., 23 (1997) 3]. Because property distributions of HDO, HAC, MW and CLOGP (used in the 'rule of 5' test) do not differ significantly between these datasets, the 'rule of 5' does not distinguish 'drugs' from 'nondrugs'. Therefore, Pareto analyses were performed to examine skewed distributions in all compound collections. Seventy percent of the 'drug-like' compounds were found between the following limits: 0 < or = HDO < or = 2, 2 < or = HAC < or = 9, 2 < or = RTB < or = 8, and 1 < or = RNG < or = 4, respectively. The number of launched drugs in MDDR having 0 < or = HDO < or = 2 is 4.8 times higher than the number of drugs having 3 < or = HDO < or = 5. Skewed distributions can be exploited to focus on the 'drug-like space': 62.68% of ACDF ('nondrug-like') compounds have 0 < or = RNG < or = 2, and RGB < or = 17, while 28.88% of ACDF compounds have 3 < or = RNG < or = 13, and 18 < or = RGB < or = 56. By contrast, 61.22% of MDDRF compounds have RNG > or = 3, and RGB > or = 18, and only 24.73% of MDDRF compounds have 0 < or = RNG < or = 2 rings, and RGB < or = 17. The probability of identifying 'drug-like' structures increases with molecular complexity.

Chemical information management in drug discovery: optimizing the computational and combinatorial chemistry interfaces.

Structure-property relationships, central to many of today's drug discovery strategies, are not straightforward to deal with when trying to predict drug efficacy, that is, the combined outcome of target affinity, pharmacodynamic behavior, pharmacokinetic properties, and metabolic fate. In this article, we discuss the handling of chemical property information in reagents-for-synthesis selection, enumeration, and virtual library construction. We describe the use of diversity assessment and/or experimental design in selection of compound-libraries-to-be-synthesized. Our overall objective was to identify good-quality drug candidates through reliable structure-activity relationship data, with the minimum number of compounds synthesized and tested. Chemical filters, property filters, scoring functions, and utilization of interactive visualization tools are discussed. The concept of chemical diversity and aspects of chemical space navigation employing a proprietary tool, Chemical Global Positioning System (ChemGPS), for mapping the drug-related chemical space are examined. Guidelines and workflow recommendations for the practicing medicinal chemist are proposed.

Toward minimalistic modeling of oral drug absorption.

Poor intestinal permeability of drugs constitutes a major bottleneck in the successful development of candidate drugs. Fast computational tools to help in designing compounds with increased probability of oral absorption are required, since both medicinal and combinatorial chemists are under pressure to consider increasing numbers of virtual and existing compounds. The QSAR paradigm for drug absorption is expressed as a function of molecular size, hydrogen-bonding capacity, and lipophilicity. A nonlinear PLS model that can be achieved with minimal computational efforts is described. The QSAR model correlates human intestinal absorption (%HIA) data, and apparent Caco-2 cell permeability data, to parameters calculated from molecular structures. Two properties were found to be relevant for absorption predictions, namely H-bonding capacity, and hydrophobic transferability. The parsimony principle was applied in several aspects: single conformers were used to compute molecular surface areas; the definitions of "polar" and "nonpolar" surfaces were done in a simplistic fashion; simple and fast 2D descriptors were used to estimate other properties; the 1 PLS component model was selected. These choices result in a minimalistic model for oral absorption. The use of both %HIA and Caco-2 permeability data was found to stabilize and improve the model. This QSAR model can serve as a simple, quantitative extension of the "rule of five" scheme (Lipinski, C.A., Lombardo, F., Dominy, B.W., and Feeney, P.J. Adv. Drug Deliv. Rev. 1997, 23, 3-25), in a manner that can prove beneficial to the drug discovery process.

The Design of Leadlike Combinatorial Libraries.

The optimization of low-potency leads into drugs is often accompanied by an increase in molecular weight (M(r)) and lipophilicity, as a consequence of affinity enhancement. Hits with affinity at µM levels discovered by screening leadlike libraries allow scope for this optimization process, as shown schematically by the distributions of M(r) for a leadlike library (1), oral drugs (2), and a typical combinatorial chemistry library (3). y=percentage with a particular molecular weight.

MTD-ADJ: a multiconformational minimal topologic difference for determining bioactive conformers using adjusted biological activities.

The active conformation is part of a conformational mixture with experimental activity Yexp, and is used in QSAR studies to extract more information regarding the ligand-receptor interaction. To reflect the relative amount (alpha) of the active conformation, we adjust Yexp: Yadj = Yexp - log alpha. We establish a quantitative structure-activity relationship (QSAR) between Yadj and 3D conformational characteristics for the acetylcholinesterase (AChE) hydrolysis rates of 25 acetic esters. The 3D-QSAR model was obtained using the adjusted multiconformational minimal steric/topologic difference (MTD-ADJ) method, optimizing the receptor map based on Yadj for each conformer. Yadj was updated during each step of the optimization process. alpha and Yadj are based on the Boltzmann distribution calculated using AMI (MOPAC 6.0) relative energies of the COSMIC 90 derived conformers. The MTD-ADJ results are: (i) the 3D-QSAR models obtained by this procedure have significant statistical parameters and are similar to the unadjusted (MTD-MC, using Yexp) models; (ii) the selected bioactive conformations are extended, occupy cavity vertices and, for the same structures, have the same MTD value; and (iii) the optimized conformational map of the neutral ligands obtained from the MTD-ADJ model fits well in the active site of the crystallographic structure of AChE (from Torpedo californica). We propose a neutral ligands binding site model for AChE. Our results show that MTD-ADJ, which can be implemented in any 3D-QSAR method, is capable of providing additional information regarding the active conformations, and can be used to gain further insight into the ligand-receptor models for which no structural data are available.

Identification of a functional water channel in cytochrome P450 enzymes.

Cytochrome P450 enzymes are monooxygenases that contain a functional heme b group linked to a conserved cysteine with a thiolate bond. In the native state, the central iron atom is hexacoordinated with a covalently bound water molecule. The exclusion of solvent molecules from the active site is essential for efficient enzymatic function. Upon substrate binding, water has to be displaced from the active site to prevent electron uncoupling that results in hydrogen peroxide or water. In contrast to typical hemoproteins, the protein surface is not directly accessible from the heme of cytochromes P450. We postulate a two-state model in which a conserved arginine, stabilizing the heme propionate in all known cytochrome P450 crystal structures, changes from the initial, stable side-chain conformation to another rotamer (metastable). In this new state, a functional water channel (aqueduct) is formed from the active site to a water cluster located on the thiolate side of the heme, close to the protein surface. This water cluster communicates with the surface in the closed state and is partly replaced by the flipping arginine side chain in the open state, allowing water molecules to exit to the surface or to reaccess the active site. This two-state model suggests the presence of an exit pathway for water between the active site and the protein surface.

Ligand-based identification of environmental estrogens.

Comparative molecular field analysis (CoMFA), a three-dimensional quantitative structure-activity relationship (3D-QSAR) paradigm, was used to examine the estrogen receptor (ER) binding affinities of a series of structurally diverse natural, synthetic, and environmental chemicals of interest. The CoMFA/3D-QSAR model is statistically robust and internally consistent, and successfully illustrates that the overall steric and electrostatic properties of structurally diverse ligands for the estrogen receptor are both necessary and sufficient to describe the binding affinity. The ability of the model to accurately predict the ER binding affinity of an external test set of molecules suggests that structure-based 3D-QSAR models may be used to supplement the process of endocrine disruptor identification through prioritization of novel compounds for bioassay. The general application of this 3D-QSAR model within a toxicological framework is, at present, limited only by the quantity and quality of biological data for relevant biomarkers of toxicity and hormonal responsiveness.

Three-dimensional model of a selective theophylline-binding RNA molecule.

A three-dimensional (3D) model for an RNA molecule that selectively binds theophylline but not caffeine is proposed. This RNA, which was found using SELEX (Jenison et al., 1994), is 10,000 times more specific for theophylline (Kn = 320 nM) than for caffeine (KD = 3.5 mM), although the two ligands are identical except for a methyl group substituted at N7 (present only in caffeine). The binding affinity for ten xanthine-based ligands was used to derive a comparative molecular field analysis model (R2 = 0.93 for three components, with cross-validated R2 of 0.73), using the SYBYL and GOLPE programs. A pharmacophoric map was generated to locate steric and electrostatic interactions between theophylline and the RNA binding site. This information was used to identify putative functional groups of the binding pocket and to generate distance constraints. On the basis of a model for the secondary structure (Jenison et al., 1994), the 3D structure of this RNA was then generated using the following method: each helical region of the RNA molecule was treated as a rigid body; single-stranded loops with specific end-to-end distances were generated. The structures of RNA-xanthine complexes were studied using a modified Monte Carlo algorithm. The detailed structure of an RNA-ligand complex model, as well as possible explanations for the theophylline selectivity are discussed.

Three-dimensional quantitative structure-activity relationships of steroid aromatase inhibitors.

Inhibition of aromatase, a cytochrome P450 that converts androgens to estrogens, is relevant in the therapeutic control of breast cancer. We investigate this inhibition using a three-dimensional quantitative structure-activity relationship (3D QSAR) method known as Comparative Molecular Field Analysis, CoMFA [Cramer III, R.D. et al., J. Am. Chem. Soc., 110 (1988) 5959]. We analyzed the data for 50 steroid inhibitors [Numazawa, M. et al., J. Med. Chem., 37 (1994) 2198, and references cited therein] assayed against androstenedione on human placental microsomes. An initial CoMFA resulted in a three-component model for log(l/Ki), with an explained variance r2 of 0.885, and a cross-validated q2 of 0.673. Chemometric studies were performed using GOLPE [Baroni, M. et al., Quant. Struct.-Act. Relatsh., 12 (1993) 9]. The CoMFA/GOLPE model is discussed in terms of robustness, predictivity, explanatory power and simplicity. After randomized exclusion of 25 or 10 compounds (repeated 25 times), the q2 for one component was 0.62 and 0.61, respectively, while r2 was 0.674. We demonstrate that the predictive r2 based on the mean activity (Ym) of the training set is misleading, while the test set Ym-based predictive r2 index gives a more accurate estimate of external predictivity. Using CoMFA, the observed differences in aromatase inhibition among C6-substituted steroids are rationalized at the atomic level. The CoMFA fields are consistent with known, potent inhibitors of aromatase, not included in the model. When positioned in the same alignment, these compounds have distinct features that overlap with the steric and electrostatic fields obtained in the CoMFA model. The presence of two hydrophobic binding pockets near the aromatase active site is discussed: a steric bulk tolerant one, common for C4, C6-alpha and C7-alpha substituents, and a smaller one at the C6-beta region.

Three-dimensional quantitative structure-activity relationship of human immunodeficiency virus (I) protease inhibitors. 2. Predictive power using limited exploration of alternate binding modes.

NewPred, a semiautomated procedure to evaluate alternate binding modes and assist three dimensional quantitative structure-activity relationship (3D-QSAR) studies in predictive power evaluation is exemplified with a series of 30 human immunodeficiency virus 1 protease (HIV PR) inhibitors. Five comparative molecular field analysis (CoMFA) models (Waller, C. L.; et al. J. Med. Chem. 1993, 36, 4152-4160) based on 59 HIV-PR inhibitors were tested. The test set included 18 compounds (set A) having a different transition state isostere (TSI), hydroxyethylurea (Getman, D. P.; et al. J. Med. Chem. 1993, 36, 288-291), to investigate the binding mode in P1' and P2'. Twelve dihyroxyethylenes (set B) (Thaisrivongs, S.; et al. J. Med. Chem. 1993, 36, 941-952) were used to investigate binding in P2 and P3 as well as in P2' and P3'. Six other compounds with known or inferred binding structure (set C) were part of the test set, but not investigated with NewPred. Each compound was aligned in accordance to predefined alignment rules for the training set prior to the inclusion in the test set (except for set C). Using NewPred, geometrically different conformers for each compound were generated and individually relaxed in the HIV-PR binding site. Energy comparisons allowed selection of lowest energy structures to be included in the test set. Only in vacuo minimized conformers derived from low-energy complexes were used to determine the predictive power of the five models (predictive r2 varied from 0.1 to 0.7 when two chemical and statistical outliers were excluded). Our models correctly predict the poor inhibitor activity of 1(S)-amino-2(R)-hydroxyindan-containing peptides (set B), which is explained and interpreted from a 3D-QSAR perspective. The use of a new, flexibility-based, semiautomated method to explore alternate binding models for 3D-QSAR models is demonstrated.