Abstract Bending of double-stranded (ds) DNA plays a crucial role in many important biological processes and is relevant for nanotechnological applications. Among all the elements that have been studied in relation to dsDNA bending, A-tracts stand out as one of the most controversial. The “ApA wedge” theory was disproved when a series of linear polynucleotides containing phased 5′-A4T4-3′ or 5′-T4A4-3′ runs were shown to be bent or straight, respectively, and crystallographic evidence revealed that A-tracts are unbent. Furthermore, some of the smallest dsDNA minicircles described to date (~ 100 bp in size) lack A-tracts and are subjected to varying levels of torsional stress. Representative DNA sequences from this experimental background were modeled in atomic detail and their dynamic behavior was simulated over hundreds of nanoseconds using the AMBER force field ParmBSC1. Subsequent analysis of the resulting trajectories allowed us to (i) unambiguously establish the location of the bends in all cases; (ii) identify the structural elements that are directly responsible for the macroscopically detected curvature; and (iii) reveal the importance not only of coherently summing the effects of the bending elements when they are in synchrony with the natural repeat of the helix (i.e. separated by an integral number of helical turns) but also when alternated with a half-integral separation of opposite effects. We conclude that the major determinant of the macroscopically observed bending is the proper grouping and phasing of the positive roll imposed by pyrimidine-purine (YR) steps and the negative or null roll characteristic of RY steps and A-tracts, respectively. This conclusion is in very good agreement with the structural parameters experimentally derived for much smaller DNA molecules either alone or as found in DNA–protein complexes. We expect that this work will pave the way for future studies on drug-induced DNA bending, DNA shape readout by transcription factors, structure of circular extrachromosomal DNA, and custom design of curved DNA origami scaffolds.

In this study quantum mechanical methods were used to predict the solvation energies of a series of drug-like molecules both in water and in octanol, in the context of the SAMPL6 n-octanol/water partition coefficient challenge. In pharmaceutical design, n-octanol/water partition coefficient, LogP, describes the drug’s hydrophobicity and membrane permeability, thus, a well-established theoretical method that rapidly determines the hydrophobicity of a drug, enables the progress of the drug design. In this study, the solvation free energies were obtained via six different methodologies (B3LYP, M06-2X and ωB97XD functionals with 6-311+G** and 6-31G* basis sets) by taking into account the environment implicitly; the methodology chosen (B3LYP/6-311+G**) was used later to evaluate ΔGsolv by using explicit water as solvent. We optimized each conformer in different solvents separately, our calculations have shown that the stability of the conformers is highly dependent on the solvent environment. We have compared implicitly and explicitly solvated systems, the interaction of one explicit water with drug-molecules at the proper location leads to the prediction of more accurate LogP values.

Abstract This work presents a quantum mechanical model for predicting octanol-water partition coefficients of small protein-kinase inhibitor fragments as part of the SAMPL6 LogP Prediction Challenge. The model calculates solvation free energy differences using the M06-2X functional with SMD implicit solvation and the def2-SVP basis set. This model was identified as dqxk4 in the SAMPL6 Challenge and was the third highest performing model in the physical methods category with 0.49 log Root Mean Squared Error (RMSE) for predicting the 11 compounds in SAMPL6 blind prediction set. We also collaboratively investigated the use of empirical models to address model deficiencies for halogenated compounds at minimal additional computational cost. A mixed model consisting of the dqxk4 physical and hdpuj empirical models found improved performance at 0.34 log RMSE on the SAMPL6 dataset. This collaborative mixed model approach shows how empirical models can be leveraged to expediently improve performance in chemical spaces that are difficult for ab initio methods to simulate.

Effective representation of a molecule is required to develop useful quantitative structure–property relationships (QSPR) for accurate prediction of chemical properties. The octanol–water partition coefficient logP, a measure of lipophilicity, is an important property for pharmacological and toxicological endpoints used in the pharmaceutical and regulatory spheres. We compare physicochemical descriptors, structural keys, and circular fingerprints in their ability to effectively represent a chemical space and characterise molecular features to correlate with lipophilicity. Exploratory landscape continuity analyses revealed that whole-molecule physicochemical descriptors could map together compounds that were similar in both molecular features and logP, indicating higher potential for use in logP QSPRs compared to the substructural approach of structural keys and circular fingerprints. Indeed, logP QSPR models parameterised by physicochemical descriptors consistently performed with the lowest error. Our best performing model was a stochastic gradient descent-optimised multilinear regression with 1438 descriptors, returning an internal benchmark RMSE of 1.03 log units. This corroborates the well-established notion that lipophilicity is an additive, whole-molecule property. We externally tested the model by participating in the 2019 SAMPL6 logP Prediction Challenge and blindly predicting for 11 protein kinase inhibitor fragment-like molecules. Our model returned an RMSE of 0.49 log units, placing eighth overall and third in the empirical methods category (submission ID ‘hdpuj’). Permutation feature importance analyses revealed that physicochemical descriptors could characterise predictive molecular features highly relevant to the kinase inhibitor fragment-like molecules.

Abstract Small molecules binding at any of the multiple regulatory sites on the molecular surface of a protein kinase may stabilize or disrupt the corresponding interaction, leading to consequent modulation of the kinase cellular activity. As such, each of these sites represents a potential drug target. Even targeting sites outside the immediate ATP site, the so-called exosites, may cause desirable biological effects through an allosteric mechanism. Targeting exosites can alleviate adverse effects and toxicity that is common when ATP-site compounds bind promiscuously to many other types of kinases. In this study we have identified, catalogued, and annotated all potentially druggable exosites on the protein kinase domains within the existing structural human kinome. We then priority-ranked these exosites by those most amenable to drug design. In order to identify pockets that are either consistent across the kinome, or unique and specific to a particular structure, we have also implemented a normalized representation of all pockets, and displayed these graphically. Finally, we have built a database and designed a web-based interface for users interested in accessing the 3-dimensional representations of these pockets. We envision this information will assist drug discovery efforts searching for untargeted binding pockets in the human kinome.

Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) offered a unique opportunity for designing and testing novel methodology for accurate docking and affinity prediction of ligands in an open and blinded manner. We participated in the beta-secretase 1 (BACE) Subchallenge which is comprised of cross-docking and redocking of 20 macrocyclic ligands to BACE and predicting binding affinity for 154 macrocyclic ligands. For this challenge, we developed machine learning models trained specifically on BACE. We developed a deep neural network (DNN) model that used a combination of both structure and ligand-based features that outperformed simpler machine learning models. According to the results released by D3R, we achieved a Spearman's rank correlation coefficient of 0.43(7) for predicting the affinity of 154 ligands. We describe the formulation of our machine learning strategy in detail. We compared the performance of DNN with linear regression, random forest, and support vector machines using ligand-based, structure-based, and combining both ligand and structure-based features. We compared different structures for our DNN and found that performance was highly dependent on fine optimization of the L2 regularization hyperparameter, alpha. We also developed a novel metric of ligand three-dimensional similarity inspired by crystallographic difference density maps to match ligands without crystal structures to similar ligands with known crystal structures. This report demonstrates that detailed parameterization, careful data training and implementation, and extensive feature analysis are necessary to obtain strong performance with more complex machine learning methods. Post hoc analysis shows that scoring functions based only on ligand features are competitive with those also using structural features. Our DNN approach tied for fifth in predicting BACE-ligand binding affinities.

Abstract We have developed an algorithm that enables simplified box orbital functions (SBO) to be obtained with optimized coefficients by fitting them to functions of many types. SBOs are a linear combination of radial functions useful in quantum chemistry calculations which can be spatially restricted (defined in \(0 \le r < r_{0}\) interval, and zero for \(r \ge r_{0}\)). The algorithm proposed makes it possible to obtain the optimal radius \(r_{0}\) and the coefficients of the SBOs of any number of terms from the functions to be fitted, but also allows the user to define a particular radius r and calculate the coefficients of the combination of terms of the SBOs. SBOs have proved to be useful in the calculation of molecular properties, and can reduce the complexity of the integral calculations, especially in huge chemical systems such as atomic clusters. These types of functions are also adequate for studying confined systems such as molecules in solution or big chemical systems such as atomic clusters. In addition, while carrying out the examples presented in this study we have tested the suitability of SBO functions to calculate molecular reactivity, showing that the basis functions provide results as good as the basis sets typically used for this kind of calculations.

In the construction of activity prediction models, the use of feature ranking methods is a useful mechanism for extracting information for ranking features in terms of their significance to develop predictive models. This paper studies the influence of feature rankers in the construction of molecular activity prediction models; for this purpose, a comparative study of fourteen rankings methods for feature selection was conducted. The activity prediction models were constructed using four well-known classifiers and a wide collection of datasets. The ranking algorithms were compared considering the performance of these classifiers using different metrics and the consistency of the ranked features.

We describe a new template-based method for docking flexible ligands such as macrocycles to proteins. It combines Monte-Carlo energy minimization on the manifold, a fast manifold search method, with BRIKARD for complex flexible ligand searching, and with the MELD accelerator of Replica-Exchange Molecular Dynamics simulations for atomistic degrees of freedom. Here we test the method in the Drug Design Data Resource blind Grand Challenge competition. This method was among the best performers in the competition, giving sub-angstrom prediction quality for the majority of the targets.

Partition coefficients describe the equilibrium partitioning of a single, defined charge state of a solute between two liquid phases in contact, typically a neutral solute. Octanol–water partition coefficients (\(K_{\rm ow}\)), or their logarithms (log P), are frequently used as a measure of lipophilicity in drug discovery. The partition coefficient is a physicochemical property that captures the thermodynamics of relative solvation between aqueous and nonpolar phases, and therefore provides an excellent test for physics-based computational models that predict properties of pharmaceutical relevance such as protein-ligand binding affinities or hydration/solvation free energies. The SAMPL6 Part II octanol–water partition coefficient prediction challenge used a subset of kinase inhibitor fragment-like compounds from the SAMPL6 \(\hbox {p}{K}_{{\rm a}}\) prediction challenge in a blind experimental benchmark. Following experimental data collection, the partition coefficient dataset was kept blinded until all predictions were collected from participating computational chemistry groups. A total of 91 submissions were received from 27 participating research groups. This paper presents the octanol–water log P dataset for this SAMPL6 Part II partition coefficient challenge, which consisted of 11 compounds (six 4-aminoquinazolines, two benzimidazole, one pyrazolo[3,4-d]pyrimidine, one pyridine, one 2-oxoquinoline substructure containing compounds) with log P values in the range of 1.95–4.09. We describe the potentiometric log P measurement protocol used to collect this dataset using a Sirius T3, discuss the limitations of this experimental approach, and share suggestions for future log P data collection efforts for the evaluation of computational methods.

Linezolid, an antibiotic of oxazolidinone family, is a translation inhibitor. The mechanism of its action that consists in preventing the binding of aminoacyl-tRNA to the A-site of the large subunit of a ribosome was embraced on the basis of the X-ray structural analysis of the linezolid complexes with vacant bacterial ribosomes. However, the known structures of the linezolid complexes with bacterial ribosomes poorly explain the linezolid selectivity in suppression of protein biosynthesis, depending on the amino acid sequence of the nascent peptide. In the present study the most probable structure of the linezolid complex with a E. coli ribosome in the A,A/P,P-state that is in line with the results of biochemical studies of linezolid action has been obtained by molecular dynamics simulation methods.

Abnormalities in the expression levels of EGFR/HER2 are found in many different types of human cancer; therefore, the design of dual inhibitors of EGFR/HER2 is a recognized anti-cancer strategy. Some lapatinib derivatives have been previously synthesized by modification at the methylsulfonylethylaminomethylfuryl group and biologically evaluated, demonstrating that the 2i compound shows potent inhibitory activity against EGFR/HER2-overexpressing cancer cells. In the present study, we explored the structural and energetic features that guide the molecular recognition of 2i using various EGFR/HER2 states. Molecular dynamics (MD) simulation with an MMPB(GB)SA approach was used to generate the inactive EGFR/HER2–ligand complexes. Our results corroborate that slight modification of lapatinib contributes to an increase in the affinity of the 2i compound for inactive EGFR/HER2 as compared with lapatinib compound, which is in accordance with experimental results. Comparison with previous results reveals that lapatinib and its derivative bind more strongly to the inactive than the intermediate active-inactive HER2 state. Principal component analysis allowed the observation that coupling of 2i to EGFR/HER2 is linked to a reduction in the conformational mobility, which may also contribute to the improvement in affinity observed for this compound as compared with lapatinib.

DNA is an essential target for the treatment of various pathologies, especially cancer. Hence targeting DNA double helix for alteration of its function has been attempted by several ways. Drug–DNA intercalation, one such biophysical process, could not be studied extensively as this requires significant deformation of the receptor DNA. Here we report thorough theoretical investigation of intercalation process in daunomycin–DNA interaction, by performing molecular dynamics simulations of the drug–DNA complexes for various DNA sequences, followed by Free-energy analysis and density functional theory (DFT) based studies to understand the binding preference. The classical energy based analyses indicate that the drug prefers to bind to TC/GA sequence over others. The DFT based energies of supra-molecular complexes are always contaminated with basis set superposition error (BSSE), which can be corrected by counterpoise method. This method is quite effective for systems containing two molecular fragments but is not appropriate for studying interaction between two base pair fragments and the drug intercalated between them. We have adopted an extension of the counterpoise method for BSSE corrected interaction energy calculation. These interaction energies, along with the energy penalty due to un-stacking of the base pairs, also indicate TC/GA sequence is the most preferred sequence for binding.

Escherichia coli (E. coli) K1 is the most common Gram-negative bacteria cause of neonatal meningitis. The penetration of E. coli through the blood–brain barrier is a key step of the meningitis pathogenesis. A host receptor protein, Caspr1, interacts with the E. coli virulence factor IbeA and thus facilitates bacterial penetration through the blood–brain barrier. Based on this result, we have now predicted the binding pattern between Caspr1 and IbeA by an integrated computational protocol. Based on the predicted model, we have identified a putative molecular binding pocket in IbeA, that directly bind with Caspr1. This evidence indicates that the IbeA (229–343aa) region might play a key role in mediating the bacteria invasion. Virtual screening with the molecular model was conducted to search for potential inhibitors from 213,279 commercially available chemical compounds. From the top 50 identified compounds, 9 demonstrated a direct binding ability to the residues within the Caspr1 binding site on IbeA. By using human brain microvascular endothelial cells (hBMEC) with E. coli strain RS218, four molecules were characterized that significantly attenuated the bacteria invasions at concentrations devoid of cell toxicity. Our study provides useful structural information for understanding the pathogenesis of neonatal meningitis, and have identified drug-like compounds that could be used to develop effective anti-meningitis agents.

Attachment of envelope glycoprotein gp120 to the host cell receptor CD4 is the first step during the human immunodeficiency virus-1 (HIV-1) entry into the host cells that makes it a promising target for drug design. To elucidate the crucial three dimensional (3D) structural features of reported HIV-1 gp120 CD4 binding inhibitors, 3D pharmacophores were generated and receptor based approach was employed to quantify these structural features. A four-partial least square factor model with good statistics and predictive ability was generated for the dataset of 100 molecules. To further ascertain the structural requirement for gp120-CD4 binding inhibition, molecular interaction studies of inhibitors with gp120 was carried out by performing molecular docking using Glide 5.6. Based on these studies, structural requirements were drawn and new molecules were designed accordingly to yield new sulphonamides derivatives. A water based green synthetic approach was adopted to obtain these compounds which were evaluated for their HIV-1 gp120 CD4 binding inhibition. The newly synthesized compounds exhibited remarkable activity (10-fold increase) when compared with the standard BMS 806. Further the stability of newly synthesized derivatives with HIV-1 gp120 was also investigated through molecular dynamics simulation studies. This provides a proof of concept for molecular modeling based design of new inhibitors for inhibition of HIV-1 gp120 CD4 interaction.

Small molecules with multi-target activity, also termed promiscuous compounds, are increasingly considered for pharmaceutical applications. The use of promiscuous chemical entities represents a departure from the compound specificity paradigm, one of the pillars of modern drug discovery. The popularity of promiscuous compounds is due to the concept of polypharmacology; another more recent drug discovery paradigm. It refers to insights that the efficacy of drugs often depends on interactions with multiple targets. Views concerning the extent to which small molecules might form well-defined interactions with multiple targets often differ, but comprehensive experimental investigations of promiscuity are currently rare. On the other hand, large volumes of active compounds and experimental measurements are becoming available and enable data-driven analyses of compound selectivity versus promiscuity. In this perspective, we discuss computational methods and data structures designed for promiscuity analysis. In addition, findings from large-scale exploration of activity profiles of inhibitors covering the human kinome are summarized. Although many kinase inhibitors are expected to be promiscuous, they are frequently found to be selective, which provides opportunities for target-directed drug discovery (rather than polypharmacology). We also discuss that machine learning yields evidence for the existence of structure–promiscuity relationships.

The D3R Grand Challenge 4 provided a brilliant opportunity to test macrocyclic docking protocols on a diverse high-quality experimental data. We participated in both pose and affinity prediction exercises. Overall, we aimed to use an automated structure-based docking pipeline built around a set of tools developed in our team. This exercise again demonstrated a crucial importance of the correct local ligand geometry for the overall success of docking. Starting from the second part of the pose prediction stage, we developed a stable pipeline for sampling macrocycle conformers. This resulted in the subangstrom average precision of our pose predictions. In the affinity prediction exercise we obtained average results. However, we could improve these when using docking poses submitted by the best predictors. Our docking tools including the Convex-PL scoring function are available at https://team.inria.fr/nano-d/software/.

The Drug Design Data Resource (D3R) Grand Challenges present an opportunity to assess, in the context of a blind predictive challenge, the accuracy and the limits of tools and methodologies designed to help guide pharmaceutical drug discovery projects. Here, we report the results of our participation in the D3R Grand Challenge 4 (GC4), which focused on predicting the binding poses and affinity ranking for compounds targeting the \(\beta\)-amyloid precursor protein (BACE-1). Our ligand similarity-based protocol using HYBRID (OpenEye Scientific Software) successfully identified poses close to the native binding mode for most of the ligands with less than 2 Å RMSD accuracy. Furthermore, we compared the performance of our HYBRID-based approach to that of AutoDock Vina and DOCK 6 and found that using a reference ligand to guide the docking process is a better strategy for pose prediction and helped HYBRID to perform better here. We also conducted end-point free energy estimates on molecules dynamics based ensembles of protein-ligand complexes using molecular mechanics combined with generalized Born surface area method (MM-GBSA). We found that the binding affinity ranking based on MM-GBSA scores have poor correlation with the experimental values. Finally, the main lessons from our participation in D3R GC4 are: (i) the generation of the macrocyclic conformers is a key step for successful pose prediction, (ii) the protonation states of the BACE-1 binding site should be treated carefully, (iii) the MM-GBSA method could not discriminate well between different predicted binding poses, and (iv) the MM-GBSA method does not perform well at predicting protein–ligand binding affinities here.

Targeting of MDM2-p53 protein–protein interaction is a current approach for the development of potent anticancer agents. The classical pharmacophore hypothesis for the design of such molecules describes the three point binding of a small molecule inhibitor to the MDM2 protein. However, this hypothesis is not confirmed when considering the activity of a number of known potent MDM2 inhibitors. Here we demonstrate the important role of the flexible N-terminal region of the MDM2 protein in the binding with small molecule compounds, which contributes to the transition from three point binding to four point binding during the development of new anticancer agents. To evaluate the contribution of the MDM2 N-terminal region to the structure–activity relationship of known MDM2 inhibitors, compounds of nutlin series, whose spatial configuration was shown to dramatically affect the target activity, were used as objects of the study. The key amino acid residues within the N-terminal region involved in the interaction with small molecule ligands were determined by means of molecular dynamics. The conformational stability of the flexible MDM2 fragment was simulated under different conditions. The effects of point mutations on the N-terminal region stability were also demonstrated.

Accurate solvation free energy ΔGsolv predictions require well parametrized force fields. In order to refit Lennard-Jones (LJ) parameters for improved ΔGsolv predictions for a variety of compound classes and chemical environments, a large number of ΔGsolv calculations is required. As the calculation of ΔGsolv is computational expensive, there is need for efficient but precise calculation methods. In this work, we focus on the computation of non-aqueous solvation free energies. We compare ΔGsolv results from highly precise reference simulations for transferring a solute from the vacuum into a condensed phase and results obtained from a thermodynamic cycle implementation. As test systems, we alter LJ parameters ε and σ of widely used GAFF atom types ca, cl, n1, oh and os in various solute/solvent combinations. We examine the degree of configurational space overlap and find an impact by hydrogen bonds and the solvent accessible surface area. We conclude that the application of a thermodynamic cycle for the parametrization of force fields targeting ΔGsolv is useful if the adaptation of LJ parameters is limited to atom types in the solute or if only ε is changed.

The IEFPCM/MST continuum solvation model is used for the blind prediction of n-octanol/water partition of a set of 11 fragment-like small molecules within the SAMPL6 Part II Partition Coefficient Challenge. The partition coefficient of the neutral species (log P) was determined using an extended parametrization of the B3LYP/6-31G(d) version of the Miertus–Scrocco–Tomasi continuum solvation model in n-octanol. Comparison with the experimental data provided for partition coefficients yielded a root-mean square error (rmse) of 0.78 (log P units), which agrees with the accuracy reported for our method (rmse = 0.80) for nitrogen-containing heterocyclic compounds. Out of the 91 sets of log P values submitted by the participants, our submission is within those with an rmse < 1 and among the four best ranked physical methods. The largest errors involve three compounds: two with the largest positive deviations (SM13 and SM08), and one with the largest negative deviations (SM15). Here we report the potentiometric determination of the log P for SM13, leading to a value of 3.62 ± 0.02, which is in better agreement with most empirical predictions than the experimental value reported in SAMPL6. In addition, further inclusion of several conformations for SM08 significantly improved our results. Inclusion of these refinements led to an overall error of 0.51 (log P units), which supports the reliability of the IEFPCM/MST model for predicting the partitioning of neutral compounds.

Imbalanced datasets, comprising of more inactive compounds relative to the active ones, are a common challenge in ligand-based model building workflows for drug discovery. This is particularly true for neglected tropical diseases since efforts to identify therapeutics for these diseases are often limited. In this report, we analyze the performance of several undersampling strategies in modeling the Dengue Virus 2 (DENV2) inhibitory activity, as well as the anti-flaviviral activities for the West Nile (WNV) and Zika (ZIKV) viruses. To this end, we build datasets comprising of 1218 (159 actives and 1059 inactives), 1044 (132 actives and 912 inactives) and 302 (75 actives and 227 inactives) molecules with known DENV2, WNV and ZIKV inhibitory activity profiles, respectively. We develop ensemble classifiers for these endpoints and compare the performance of the different undersampling algorithms on external sets. It is observed that data pruning algorithms yield superior performance relative to data selection algorithms. The best overall performance is provided by the one-sided selection algorithm with test set balanced accuracy (BACC) values of 0.84, 0.74 and 0.77 for the DENV2, WNV and ZIKV inhibitory activities, respectively. For the model building, we use the recently proposed GT-STAF information indices, and compare the predictivity of 3 molecular fragmentation approaches: connected subgraphs, substructure and alogp atom types, which are observed to show comparable performance. On the other hand, a combination of indices based on these fragmentation strategies enhances the predictivity of the built ensembles. The built models could be useful for screening new molecules with possible DENV, WNV and ZIKV inhibitory activities. ADMET modelers are encouraged to adopt undersampling algorithms in their workflows when dealing with imbalanced datasets.

Influenza epidemics are responsible for an average of 3–5 millions of severe cases and up to 500,000 deaths around the world. One of flu pandemic types is influenza A(H1N1)pdm09 virus (pdm09H1N1). Oseltamivir is the antiviral drug used to treat influenza targeting at neuraminidase (NA) located on the viral surface. Influenza virus undergoes high mutation rates and leads to drug resistance, and thus the development of more efficient drugs is required. In the present study, all-atom molecular dynamics simulations were applied to understand the oseltamivir resistance caused by the single E119D and double E119D/H274Y mutations on NA. The obtained results in terms of binding free energy and intermolecular interactions in the ligand–protein interface showed that the oseltamivir could not be well accommodated in the binding pocket of both NA mutants and the 150-loop moves out from oseltamivir as an “open” state. A greater number of water molecules accessible to the binding pocket could disrupt the oseltamivir binding with NA target as seen be high mobility of oseltamivir at the active site. Additionally, our finding could guide to the design and development of novel NA inhibitor drugs.

Within the framework of the 6th physical property blind challenge (SAMPL6) the authors have participated in predicting the octanol–water partition coefficients (logP) for several small drug like molecules. Those logP values where experimentally known by the organizers but only revealed after the submissions of the predictions. Two different sets of predictions were submitted by the authors, both based on the COSMOtherm implementation of COSMO-RS theory. COSMOtherm predictions using the FINE parametrization level (hmz0n) obtained the highest accuracy among all submissions as measured by the root mean squared error. COSMOquick predictions using a fast algorithm to estimate σ-profiles and an a posterio machine learning correction on top of the COSMOtherm results (3vqbi) scored 3rd out of 91 submissions. Both results underline the high quality of COSMO-RS derived molecular free energies in solution.

Application of a small radius of the hydrogen atom in molecular-mechanical models of hydrogen bonding improves the estimate of the solvation free energy of organic substances. At the same time, the density and evaporation heat of the bulk water vary slightly and are close to experimental values. Blind testing drug candidates in the SAMPL6 simulation competition showed that using the same Lennard-Jones hydrogen parameters for the hydroxyl, hydroxycarbonyl, amino, and amide groups is enough to predict log P octanol–water with MAE 0.67 and RMSD 0.75.

We investigated the Gi-coupled A3 adenosine receptor (A3AR) activation mechanism by running 7.2 µs of molecular dynamics (MD) simulations. Based on homology to G protein-coupled receptor (GPCR) structures, three constitutively active mutant (CAM) and the wild-type (WT) A3ARs in the apo form were modeled. Conformational signatures associated with three different receptor states (inactive R, active R*, and bound to Gi protein mimic) were predicted by analyzing and comparing the CAMs with WT receptor and by considering site-directed mutagenesis data. Detected signatures that were correlated with receptor state included: Persistent salt-bridges involving key charged residues for activation (including a novel, putative ionic lock), rotameric state of conserved W6.48, and Na+ ions and water molecules present. Active-coupled state signatures similar to the X-ray structures of β2 adrenergic receptor-Gs protein and A2AAR-mini-Gs and the recently solved cryo-EM A1AR–Gi complexes were found. Our MD analysis suggests that constitutive activation might arise from the D1073.49–R1083.50 ionic lock destabilization in R and the D1073.49–R1113.53 ionic lock stabilization in R* that presumably lowers the energy barrier associated with an R to R* transition. This study provides new opportunities to understand the underlying interactions of different receptor states of other Gi protein-coupled GPCRs.

The Glucagon-like peptide 1 receptor (GLP-1R) is a well-established target for the treatment of type 2 diabetes and GLP-1R agonist-based therapies represent an effective approach which results in several GLP-1 analog drugs. However, the development of nonpeptidic agonist drugs targeting GLP-1R remains unsuccessful. A promising strategy aims to develop orally bioavailable, small-molecule positive allosteric modulators of GLP1-1R. Taking advantage of the recently reported cryo-EM structure of GLP-1R at its active state, we have performed structure-based screening studies which include potential allosteric binding site prediction and in silico screening of drug-like compounds, and conducted in vitro testing and site-specific mutagenesis studies. One compound with low molecular weight was confirmed as a positive allosteric modulator of GLP-1R as it enhances GLP-1′s affinity and efficacy to human GLP-1R in a dose dependent manner. This compound also stimulates insulin secretion synergistically with GLP-1. With the molecular weight of 399, this compound represents one of the smallest known GLP-1R PAMs, and demonstrates other favorable drug-like properties. Site-specific mutagenesis studies confirmed that the binding site of this compound partially overlaps with that of a known antagonist in the transmembrane domain. These results demonstrate that structure-based approach is useful for discovering nonpeptidic allosteric modulators of GLP-1R and the compound reported here is valuable for further drug development.

Development of novel in silico methods for questing novel PgP inhibitors is crucial for the reversal of multi-drug resistance in cancer therapy. Here, we report machine learning based binary classification schemes to identify the PgP inhibitors from non-inhibitors using molecular solvation theory with excellent accuracy and precision. The excess chemical potential and partial molar volume in various solvents are calculated for PgP± (PgP inhibitors and non-inhibitors) compounds with the statistical–mechanical based three-dimensional reference interaction site model with the Kovalenko–Hirata closure approximation (3D-RISM-KH molecular theory of solvation). The statistical importance analysis of descriptors identified the 3D-RISM-KH based descriptors as top molecular descriptors for classification. Among the constructed classification models, the support vector machine predicted the test set of Pgp± compounds with highest accuracy and precision of ~ 97% for test set. The validation of models confirms the robustness of state-of-the-art molecular solvation theory based descriptors in identification of the Pgp± compounds.

The in silico prediction of partition coefficients is an important task in computer-aided drug discovery. In particular the octanol–water partition coefficient is used as surrogate for lipophilicity. Various computational approaches have been proposed, ranging from simple group-contribution techniques based on the 2D topology of a molecule to rigorous methods based molecular dynamics (MD) or quantum chemistry. In order to balance accuracy and computational cost, we recently developed the MD fingerprints (MDFPs), where the information in MD simulations is encoded in a floating-point vector, which can be used as input for machine learning (ML). The MDFP-ML approach was shown to perform similarly to rigorous methods while being substantially more efficient. Here, we present the application of MDFP-ML for the prediction of octanol–water partition coefficients in the SAMPL6 blind challenge. The underlying computational pipeline is made freely available in form of the MDFPtools package.

We present the performances of our mathematical deep learning (MathDL) models for D3R Grand Challenge 4 (GC4). This challenge involves pose prediction, affinity ranking, and free energy estimation for beta secretase 1 (BACE) as well as affinity ranking and free energy estimation for Cathepsin S (CatS). We have developed advanced mathematics, namely differential geometry, algebraic graph, and/or algebraic topology, to accurately and efficiently encode high dimensional physical/chemical interactions into scalable low-dimensional rotational and translational invariant representations. These representations are integrated with deep learning models, such as generative adversarial networks (GAN) and convolutional neural networks (CNN) for pose prediction and energy evaluation, respectively. Overall, our MathDL models achieved the top place in pose prediction for BACE ligands in Stage 1a. Moreover, our submissions obtained the highest Spearman correlation coefficient on the affinity ranking of 460 CatS compounds, and the smallest centered root mean square error on the free energy set of 39 CatS molecules. It is worthy to mention that our method on docking pose predictions has significantly improved from our previous ones.

Cathepsin S (CatS), a member of cysteine cathepsin proteases, has been well studied due to its significant role in many pathological processes, including arthritis, cancer and cardiovascular diseases. CatS inhibitors have been included in D3R-GC3 for both docking pose prediction and affinity ranking, and in D3R-GC4 for binding affinity ranking. The difficulties posed by CatS inhibitors in D3R mainly come from three aspects: large size, high flexibility and similar chemical structures. We have participated in GC4; our best submitted model, which employs a similarity-based alignment docking and Vina scoring protocol, yielded Kendall’s τ of 0.23 for 459 binders in GC4. In our further explorations with machine learning, by curating a CatS specific training set, adopting a similarity-based constrained docking method as well as an arm-based fragmentation strategy which can describe large inhibitors in a locality-sensitive fashion, our best structure-based ranking protocol can achieve Kendall’s τ of 0.52 for all binders in GC4. In this exploration process, we have demonstrated the importance of training data, docking approaches and fragmentation strategies in inhibitor-ranking protocol development with machine learning.

The computational prediction of ligand-biopolymer affinities is a crucial endeavor in modern drug discovery and one that still poses major challenges. The choice of the appropriate computational method often reveals itself as a trade-off between accuracy and speed, with mathematical devices referred to as scoring functions being the fastest. Among the many shortcomings of scoring functions there is the lack of universal applicability to every molecular system. This is so largely due to their reliance on atom type perception and/or parametrization. This article proposes the use of nonparametric Model of Effective Radii of Atoms descriptors that can be readily computed for the entire Periodic Table and demonstrate that, in combination with machine learning algorithms, they can yield competitive performances and chemically meaningful insights.

We report the performance of our newly introduced Ensemble Docking with Enhanced sampling of pocket Shape (EDES) protocol coupled to a template-based algorithm to generate near-native ligand conformations in the 2019 iteration of the Grand Challenge (GC4) organized by the D3R consortium. Using either AutoDock4.2 or HADDOCK2.2 docking programs (each software in two variants of the protocol) our method generated native-like poses among the top 5 submitted for evaluation for most of the 20 targets with similar performances. The protein selected for GC4 was the human beta-site amyloid precursor protein cleaving enzyme 1 (BACE-1), a transmembrane aspartic-acid protease. We identified at least one pose whose heavy-atoms RMSD was less than 2.5 Å from the native conformation for 16 (80%) and 17 (85%) of the 20 targets using AutoDock and HADDOCK, respectively. Dissecting the possible sources of errors revealed that: (i) our EDES protocol (with minor modifications) was able to sample sub-ångstrom conformations for all 20 protein targets, reproducing the correct conformation of the binding site within ~ 1 Å RMSD; (ii) as already shown by some of us in GC3, even in the presence of near-native protein structures, a proper selection of ligand conformers is crucial for the success of ensemble-docking calculations. Importantly, our approach performed best among the protocols exploiting only structural information of the apo protein to generate conformations of the receptor for ensemble-docking calculations.

The original version of this article unfortunately contained some mistakes in the references.

In this paper we describe our approaches to predict the binding mode of twenty BACE1 ligands as part of Grand Challenge 4 (GC4), organized by the Drug Design Data Resource. Calculations for all submissions (except for one, which used AutoDock4.2) were performed using AutoDock-GPU, the new GPU-accelerated version of AutoDock4 implemented in OpenCL, which features a gradient-based local search. The pose prediction challenge was organized in two stages. In Stage 1a, the protein conformations associated with each of the ligands were undisclosed, so we docked each ligand to a set of eleven receptor conformations, chosen to maximize the diversity of binding pocket topography. Protein conformations were made available in Stage 1b, making it a re-docking task. For all calculations, macrocyclic conformations were sampled on the fly during docking, taking the target structure into account. To leverage information from existing structures containing BACE1 bound to ligands available in the PDB, we tested biased docking and pose filter protocols to facilitate poses resembling those experimentally determined. Both pose filters and biased docking resulted in more accurate docked poses, enabling us to predict for both Stages 1a and 1b ligand poses within 2 Å RMSD from the crystallographic pose. Nevertheless, many of the ligands could be correctly docked without using existing structural information, demonstrating the usefulness of physics-based scoring functions, such as the one used in AutoDock4, for structure based drug design.

Molecular docking has been successfully used in computer-aided molecular design projects for the identification of ligand poses within protein binding sites. However, relying on docking scores to rank different ligands with respect to their experimental affinities might not be sufficient. It is believed that the binding scores calculated using molecular mechanics combined with the Poisson–Boltzman surface area (MM-PBSA) or generalized Born surface area (MM-GBSA) can predict binding affinities more accurately. In this perspective, we decided to take part in Stage 2 of the Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) to compare the performance of a quick scoring function, AutoDock4, to that of MM-GBSA in predicting the binding affinities of a set of \(\beta\)-Amyloid Cleaving Enzyme 1 (BACE-1) ligands. Our results show that re-scoring docking poses using MM-GBSA did not improve the correlation with experimental affinities. We further did a retrospective analysis of the results and found that our MM-GBSA protocol is sensitive to details in the protein-ligand system: (i) neutral ligands are more adapted to MM-GBSA calculations than charged ligands, (ii) predicted binding affinities depend on the initial conformation of the BACE-1 receptor, (iii) protonating the aspartyl dyad of BACE-1 correctly results in more accurate binding affinity predictions.

Mycobacterium tuberculosis infection remains a major cause of global morbidity and mortality due to the increase of antibiotics resistance. Dual/multi-target drug discovery is a promising approach to overcome bacterial resistance. In this study, we built ligand-based pharmacophore models and performed pharmacophore screening in order to identify hit compounds targeting simultaneously two enzymes—M. tuberculosis leucyl-tRNA synthetase (LeuRS) and methionyl-tRNA synthetase (MetRS). In vitro aminoacylation assay revealed five compounds from different chemical classes inhibiting both enzymes. Among them the most active compound—3-(3-chloro-4-methoxy-phenyl)-5-[3-(4-fluoro-phenyl)-[1,2,4]oxadiazol-5-yl]-3H-[1,2,3]triazol-4-ylamine (1) inhibits mycobacterial LeuRS and MetRS with IC50 values of 13 µM and 13.8 µM, respectively. Molecular modeling study indicated that compound 1 has similar binding mode with the active sites of both aminoacyl-tRNA synthetases and can be valuable compound for further chemical optimization in order to find promising antituberculosis agents.

Misfolded Cu/Zn superoxide dismutase enzyme (SOD1) shows prion-like propagation in neuronal cells leading to neurotoxic aggregates that are implicated in amyotrophic lateral sclerosis (ALS). Tryptophan-32 (W32) in SOD1 is part of a potential site for templated conversion of wild type SOD1. This W32 binding site is located on a convex, solvent exposed surface of the SOD1 suggesting that hydration effects can play an important role in ligand recognition and binding. A recent X-ray crystal structure has revealed that 5-Fluorouridine (5-FUrd) binds at the W32 binding site and can act as a pharmacophore scaffold for the development of anti-ALS drugs. In this study, a new protocol is developed to account for structural (non-displaceable) water molecules in docking simulations and successfully applied to predict the correct docked conformation binding modes of 5-FUrd at the W32 binding site. The docked configuration is within 0.58 Å (RMSD) of the observed configuration. The docking protocol involved calculating a hydration structure around SOD1 using molecular theory of solvation (3D-RISM-KH, 3D-Reference Interaction Site Model-Kovalenko-Hirata) whereby, non-displaceable water molecules are identified for docking simulations. This protocol was also used to analyze the hydrated structure of the W32 binding site and to explain the role of solvation in ligand recognition and binding to SOD1. Structural water molecules mediate hydrogen bonds between 5-FUrd and the receptor, and create an environment favoring optimal placement of 5-FUrd in the W32 binding site.

Using the D3R Grand Challenge 4 dataset containing Beta-secretase 1 (BACE) and Cathepsin S (CatS) inhibitors, we have evaluated the performance of our in-house docking workflow that involves in the first step the selection of the most suitable docking software for the system of interest based on structural and functional information available in public databases, followed by the docking of the dataset to predict the binding modes and ranking of ligands. The macrocyclic nature of the BACE ligands brought additional challenges, which were dealt with by a careful preparation of the three-dimensional input structures for ligands. This provided top-performing predictions for BACE, in contrast with CatS, where the predictions in the absence of guiding constraints provided poor results. These results highlight the importance of previous structural knowledge that is needed for correct predictions on some challenging targets. After the end of the challenge, we also carried out free energy calculations (i.e. in a non-blinded manner) for CatS using the pmx software and several force fields (AMBER, Charmm). Using knowledge-based starting pose construction allowed reaching remarkable accuracy for the CatS free energy estimates. Interestingly, we show that the use of a consensus result, by averaging the results from different force fields, increases the prediction accuracy.

We present a Focused Library Generator that is able to create from scratch new molecules with desired properties. After training the Generator on the ChEMBL database, transfer learning was used to switch the generator to producing new Mdmx inhibitors that are a promising class of anticancer drugs. Lilly medicinal chemistry filters, molecular docking, and a QSAR IC50 model were used to refine the output of the Generator. Pharmacophore screening and molecular dynamics (MD) simulations were then used to further select putative ligands. Finally, we identified five promising hits with equivalent or even better predicted binding free energies and IC50 values than known Mdmx inhibitors. The source code of the project is available on https://github.com/bigchem/online-chem.

Antagonism of CCR9 is a promising mechanism for treatment of inflammatory bowel disease, including ulcerative colitis and Crohn's disease. There is limited experimental data on CCR9 and its ligands, complicating efforts to identify new small molecule antagonists. We present here results of a successful virtual screening and rational hit-to-lead campaign that led to the discovery and initial optimization of novel CCR9 antagonists. This work uses a novel data fusion strategy to integrate the output of multiple computational tools, such as 2D similarity search, shape similarity, pharmacophore searching, and molecular docking, as well as the identification and incorporation of privileged chemokine fragments. The application of various ranking strategies, which combined consensus and parallel selection methods to achieve a balance of enrichment and novelty, resulted in 198 virtual screening hits in total, with an overall hit rate of 18%. Several hits were developed into early leads through targeted synthesis and purchase of analogs.

We have calculated relative binding affinities for eight tetrafluorophenyl-triazole-thiogalactoside inhibitors of galectin-3 with the alchemical free-energy perturbation approach. We obtain a mean absolute deviation from experimental estimates of only 2-3 kJ/mol and a correlation coefficient (R2) of 0.5-0.8 for seven relative affinities spanning a range of up to 11 kJ/mol. We also studied the effect of using different methods to calculate the charges of the inhibitor and different sizes of the perturbed group (the atoms that are described by soft-core potentials and are allowed to have differing coordinates). However, the various approaches gave rather similar results and it is not possible to point out one approach as consistently and significantly better than the others. Instead, we suggest that such small and reasonable variations in the computational method can be used to check how stable the calculated results are and to obtain a more accurate estimate of the uncertainty than if performing only one calculation with a single computational setup.

Colchicine site inhibitors are microtubule destabilizers having promising role in cancer therapeutics. In the current study, four such indanone derivatives (t1, t9, t14 and t17) with 3,4,5-trimethoxyphenyl fragment (ring A) and showing significant microtubule destabilization property have been explored. The interaction mechanism and conformational modes triggered by binding of these indanone derivatives and combretastatin at colchicine binding site (CBS) of αβ-tubulin dimer were studied using molecular dynamics (MD) simulation, principle component analysis and free energy landscape analysis. In the MD results, t1 showed binding similar to colchicine interacting in the deep hydrophobic core at the CBS. While t9, t14 and t17 showed binding conformation similar to combretastatin, with ring A superficially binding at the CBS. Results demonstrated that ring A played a vital role in binding via hydrophobic interactions and got anchored between the S8 and S9 sheets, H8 helix and T7 loop at the CBS. Conformational modes study revealed that twisting and bending conformational motions (as found in the apo system) were nearly absent in the ligand bound systems. Absence of twisting motion might causes loss of lateral contacts in microtubule, thus promoting microtubule destabilization. This study provides detailed account of microtubule destabilization mechanism by indanone ligands and combretastatin, and would be helpful for designing microtubule destabilizers with higher activity.

Acetyl-coenzyme A carboxylases (ACCs) is the first committed enzyme of fatty acid synthesis pathway. The inhibition of ACC is thought to be beneficial not only for diseases related to metabolism, such as type-2 diabetes, but also for infectious disease like bacterial infection disease. Soraphen A, a potent allosteric inhibitor of BC domain of yeast ACC, exhibit lower binding affinities to several yeast ACC mutants and the corresponding drug resistance mechanisms are still unknown. We report here a theoretical study of binding of soraphen A to wild type and yeast ACC mutants (including F510I, N485G, I69E, E477R, and K73R) via molecular dynamic simulation and molecular mechanics/generalized Born surface area free energy calculations methods. The calculated binding free energies of soraphen A to yeast ACC mutants are weaker than to wild type, which is highly consistent with the experimental results. The mutant F510I weakens the binding affinity of soraphen A to yeast ACC mainly by decreasing the van der Waals contributions, while the weaker binding affinities of Soraphen A to other yeast ACC mutants including N485G, I69E, E477R, and K73R are largely attributed to the decreased net electrostatic (ΔEele + ΔGGB) interactions. Our simulation results could provide important insights for the development of more potent ACC inhibitors.

In the development of novel pharmaceuticals, the knowledge of how many, and which, Cytochrome P450 isoforms are involved in the phase I metabolism of a compound is important. Potential problems can arise if a compound is metabolised predominantly by a single isoform in terms of drug-drug interactions or genetic polymorphisms that would lead to variations in exposure in the general population. Combined with models of regioselectivities of metabolism by each isoform, such a model would also aid in the prediction of the metabolites likely to be formed by P450-mediated metabolism. We describe the generation of a multi-class random forest model to predict which, out of a list of the seven leading Cytochrome P450 isoforms, would be the major metabolising isoforms for a novel compound. The model has a 76% success rate with a top-1 criterion and an 88% success rate for a top-2 criterion and shows significant enrichment over randomised models.

The quantitative structure-activity relationships method was popularized by Hansch and Fujita over 50 years ago. The usefulness of the method for drug design and development has been shown in the intervening years. As it was developed initially to elucidate which molecular properties modulated the relative potency of putative agrochemicals, and at a time when computing resources were scarce, there is much scope for applying modern mathematical methods to improve the QSAR method and to extending the general concept to the discovery and optimization of bioactive molecules and materials more broadly. I describe research over the past two decades where we have rebuilt the unit operations of the QSAR method using improved mathematical techniques, and have applied this valuable platform technology to new important areas of research and industry such as nanoscience, omics technologies, advanced materials, and regenerative medicine. This paper was presented as the 2017 ACS Herman Skolnik lecture.

Understanding how proteins encode ligand specificity is fascinating and similar in importance to deciphering the genetic code. For protein-ligand recognition, the combination of an almost infinite variety of interfacial shapes and patterns of chemical groups makes the problem especially challenging. Here we analyze data across non-homologous proteins in complex with small biological ligands to address observations made in our inhibitor discovery projects: that proteins favor donating H-bonds to ligands and avoid using groups with both H-bond donor and acceptor capacity. The resulting clear and significant chemical group matching preferences elucidate the code for protein-native ligand binding, similar to the dominant patterns found in nucleic acid base-pairing. On average, 90% of the keto and carboxylate oxygens occurring in the biological ligands formed direct H-bonds to the protein. A two-fold preference was found for protein atoms to act as H-bond donors and ligand atoms to act as acceptors, and 76% of all intermolecular H-bonds involved an amine donor. Together, the tight chemical and geometric constraints associated with satisfying donor groups generate a hydrogen-bonding lock that can be matched only by ligands bearing the right acceptor-rich key. Measuring an index of H-bond preference based on the observed chemical trends proved sufficient to predict other protein-ligand complexes and can be used to guide molecular design. The resulting Hbind and Protein Recognition Index software packages are being made available for rigorously defining intermolecular H-bonds and measuring the extent to which H-bonding patterns in a given complex match the preference key.

David Weininger's career, accomplishments, genius, and friendship are warmly remembered by several of his colleagues, friends, and admirers.

Antagonists of the nociceptin receptor (NOP) are raising interest for their possible clinical use as antidepressant drugs. Recently, the structure of NOP in complex with some piperidine-based antagonists has been revealed by X-ray crystallography. In this study, a multi-flexible docking (MF-docking) procedure, i.e. docking to multiple receptor conformations extracted by preliminary molecular dynamics trajectories, together with hybrid quantum mechanics/molecular mechanics (QM/MM) simulations have been carried out to provide the binding mode of two novel NOP antagonists, one of them selective (BTRX-246040, formerly named LY-2940094) and one non selective (AT-076), i.e. able to inactivate NOP as well as the classical µ- k- and δ-opioid receptors (MOP KOP and DOP). According to our results, the pivotal role of residue D1303,32 (upper indexes are Ballesteros-Weinstein notations) is analogous to that enlighten by the already known X-ray structures of opioid receptors: binding of the molecules are predicted to require a slight readjustment of the hydrophobic pocket (residues Y1313,33, M1343,36, I2195,43, Q2806,52 and V2836,55) in the orthosteric site of NOP, accommodating either the pyridine-pyrazole (BTRX-246040) or the isoquinoline (AT-076) moiety of the ligand, in turn allowing the protonated piperidine nitrogen to maximize interaction (salt-bridge) with residue D1303,32 of the NOP, and the aromatic head to be sandwiched in optimal π-stacking between Y1313,33 and M1343,36. The QM/MM optimization after the MF-docking procedure has provided the more likely conformations for the binding to the NOP receptor of BTRX-246040 and AT-076, based on different pharmacophores and exhibiting different selectivity profiles. While the high selectivity for NOP of BTRX-246040 can be explained by interactions with NOP specific residues, the lack of selectivity of AT-076 could be associated to its ability to penetrate into the deep hydrophobic pocket of NOP, while retaining a conformation very similar to the one assumed by the antagonist JDTic into the K-opioid receptor. The proposed binding geometries fit better the binding pocket environment providing clues for experimental studies aimed to design selective or multifunctional opioid drugs.

Drug-target networks have aided in many target prediction studies aiming at drug repurposing or the analysis of side effects. Conventional drug-target networks are bipartite. They contain two different types of nodes representing drugs and targets, respectively, and edges indicating pairwise drug-target interactions. In this work, we introduce a tripartite network consisting of drugs, other bioactive compounds, and targets from different sources. On the basis of analog relationships captured in the network and so-called neighbor targets of drugs, new drug targets can be inferred. The tripartite network was found to have a stable structure and simulated network growth was accompanied by a steady increase in assortativity, reflecting increasing correlation between degrees of connected nodes leading to even network connectivity. Local drug environments in the tripartite network typically contained neighbor targets and revealed interesting drug-compound-target relationships for further analysis. Candidate targets were prioritized. The tripartite network design extends standard drug-target networks and provides additional opportunities for drug target prediction.

Optimization of fragment size D-amino acid oxidase (DAAO) inhibitors was investigated using a combination of computational and experimental methods. Retrospective free energy perturbation (FEP) calculations were performed for benzo[d]isoxazole derivatives, a series of known inhibitors with two potential binding modes derived from X-ray structures of other DAAO inhibitors. The good agreement between experimental and computed binding free energies in only one of the hypothesized binding modes strongly support this bioactive conformation. Then, a series of 1-H-indazol-3-ol derivatives formerly not described as DAAO inhibitors was investigated. Binding geometries could be reliably identified by structural similarity to benzo[d]isoxazole and other well characterized series and FEP calculations were performed for several tautomers of the deprotonated and protonated compounds since all these forms are potentially present owing to the experimental pKa values of representative compounds in the series. Deprotonated compounds are proposed to be the most important bound species owing to the significantly better agreement between their calculated and measured affinities compared to the protonated forms. FEP calculations were also used for the prediction of the affinities of compounds not previously tested as DAAO inhibitors and for a comparative structure-activity relationship study of the benzo[d]isoxazole and indazole series. Selected indazole derivatives were synthesized and their measured binding affinity towards DAAO was in good agreement with FEP predictions.

PI3Kα is a promising drug target for cancer chemotherapy. In this paper, we report a strategy of combing ligand-based and structure-based virtual screening to identify new PI3Kα inhibitors. First, naïve Bayesian (NB) learning models and a 3D-QSAR pharmacophore model were built based upon known PI3Kα inhibitors. Then, the SPECS library was screened by the best NB model. This resulted in virtual hits, which were validated by matching the structures against the pharmacophore models. The pharmacophore matched hits were then docked into PI3Kα crystal structures to form ligand-receptor complexes, which are further validated by the Glide-XP program to result in structural validated hits. The structural validated hits were examined by PI3Kα inhibitory assay. With this screening protocol, ten PI3Kα inhibitors with new scaffolds were discovered with IC50 values ranging 0.44-31.25 μM. The binding affinities for the most active compounds 33 and 74 were estimated through molecular dynamics simulations and MM-PBSA analyses.

Quantitative structure-activity relationship (QSAR) is an effective computational technique for drug design that relates the chemical structures of compounds to their biological activities. Feature selection is an important step in QSAR based drug design to select the most relevant descriptors. One of the most popular feature selection methods for classification problems is Fisher score which aim is to minimize the within-class distance and maximize the between-class distance. In this study, the properties of Fisher criterion were extended for QSAR models to define the new distance metrics based on the continuous activity values of compounds with known activities. Then, a semi-supervised feature selection method was proposed based on the combination of Fisher and Laplacian criteria which exploits both compounds with known and unknown activities to select the relevant descriptors. To demonstrate the efficiency of the proposed semi-supervised feature selection method in selecting the relevant descriptors, we applied the method and other feature selection methods on three QSAR data sets such as serine/threonine-protein kinase PLK3 inhibitors, ROCK inhibitors and phenol compounds. The results demonstrated that the QSAR models built on the selected descriptors by the proposed semi-supervised method have better performance than other models. This indicates the efficiency of the proposed method in selecting the relevant descriptors using the compounds with known and unknown activities. The results of this study showed that the compounds with known and unknown activities can be helpful to improve the performance of the combined Fisher and Laplacian based feature selection methods.

In virtual drug screening, the chemical diversity of hits is an important factor, along with their predicted activity. Moreover, interim results are of interest for directing the further research, and their diversity is also desirable. In this paper, we consider a problem of obtaining a diverse set of virtual screening hits in a short time. To this end, we propose a mathematical model of task scheduling for virtual drug screening in high-performance computational systems as a congestion game between computational nodes to find the equilibrium solutions for best balancing the number of interim hits with their chemical diversity. The model considers the heterogeneous environment with workload uncertainty, processing time uncertainty, and limited knowledge about the input dataset structure. We perform computational experiments and evaluate the performance of the developed approach considering organic molecules database GDB-9. The used set of molecules is rich enough to demonstrate the feasibility and practicability of proposed solutions. We compare the algorithm with two known heuristics used in practice and observe that game-based scheduling outperforms them by the hit discovery rate and chemical diversity at earlier steps. Based on these results, we use a social utility metric for assessing the efficiency of our equilibrium solutions and show that they reach greatest values.

Use of quantum mechanical/molecular mechanical (QM/MM) methods in binding free energy calculations, particularly in the SAMPL challenge, often fail to achieve improvement over standard additive (MM) force fields. Frequently, the implementation is through use of reference potentials, or the so-called "indirect approach", and inherently relies on sufficient overlap existing between MM and QM/MM configurational spaces. This overlap is generally poor, particularly for the use of free energy perturbation to perform the MM to QM/MM free energy correction at the end states of interest (e.g., bound and unbound states). However, by utilizing MM parameters that best reproduce forces obtained at the desired QM level of theory, it is possible to lessen the configurational disparity between MM and QM/MM. To this end, we sought to use force matching to generate MM parameters for the SAMPL6 CB[8] host-guest binding challenge, classically compute binding free energies, and apply energetic end state corrections to obtain QM/MM binding free energy differences. For the standard set of 11 molecules and the bonus set (including three additional challenge molecules), error statistics, such as the root mean square deviation (RMSE) were moderately poor (5.5 and 5.4 kcal/mol). Correlation statistics, however, were in the top two for both standard and bonus set submissions ([Formula: see text] of 0.42 and 0.26, [Formula: see text] of 0.64 and 0.47 respectively). High RMSE and moderate correlation strongly indicated the presence of systematic error. Identifiable issues were ameliorated for two of the guest molecules, resulting in a reduction of error and pointing to strong prospects for the future use of this methodology.

Accurately predicting the binding affinities of small organic molecules to biological macromolecules can greatly accelerate drug discovery by reducing the number of compounds that must be synthesized to realize desired potency and selectivity goals. Unfortunately, the process of assessing the accuracy of current computational approaches to affinity prediction against binding data to biological macromolecules is frustrated by several challenges, such as slow conformational dynamics, multiple titratable groups, and the lack of high-quality blinded datasets. Over the last several SAMPL blind challenge exercises, host-guest systems have emerged as a practical and effective way to circumvent these challenges in assessing the predictive performance of current-generation quantitative modeling tools, while still providing systems capable of possessing tight binding affinities. Here, we present an overview of the SAMPL6 host-guest binding affinity prediction challenge, which featured three supramolecular hosts: octa-acid (OA), the closely related tetra-endo-methyl-octa-acid (TEMOA), and cucurbit[8]uril (CB8), along with 21 small organic guest molecules. A total of 119 entries were received from ten participating groups employing a variety of methods that spanned from electronic structure and movable type calculations in implicit solvent to alchemical and potential of mean force strategies using empirical force fields with explicit solvent models. While empirical models tended to obtain better performance than first-principle methods, it was not possible to identify a single approach that consistently provided superior results across all host-guest systems and statistical metrics. Moreover, the accuracy of the methodologies generally displayed a substantial dependence on the system considered, emphasizing the need for host diversity in blind evaluations. Several entries exploited previous experimental measurements of similar host-guest systems in an effort to improve their physical-based predictions via some manner of rudimentary machine learning; while this strategy succeeded in reducing systematic errors, it did not correspond to an improvement in statistical correlation. Comparison to previous rounds of the host-guest binding free energy challenge highlights an overall improvement in the correlation obtained by the affinity predictions for OA and TEMOA systems, but a surprising lack of improvement regarding root mean square error over the past several challenge rounds. The data suggests that further refinement of force field parameters, as well as improved treatment of chemical effects (e.g., buffer salt conditions, protonation states), may be required to further enhance predictive accuracy.

Methylphenidate (MP) binds to the cocaine binding site on the dopamine transporter and inhibits reuptake of dopamine, but does not appear to have the same abuse potential as cocaine. This study, part of a comprehensive effort to identify a drug treatment for cocaine abuse, investigates the effect of choice of calculation technique and of solvent model on the conformational potential energy surface (PES) of MP and a rigid methylphenidate (RMP) analogue which exhibits the same dopamine transporter binding affinity as MP. Conformational analysis was carried out by the AM1 and AM1/SM5.4 semiempirical molecular orbital methods, a molecular mechanics method (Tripos force field with the dielectric set equal to that of vacuum or water) and the HF/6-31G* molecular orbital method in vacuum phase. Although all three methods differ somewhat in the local details of the PES, the general trends are the same for neutral and protonated MP. In vacuum phase, protonation has a distinctive effect in decreasing the regions of space available to the local conformational minima. Solvent has little effect on the PES of the neutral molecule and tends to stabilize the protonated species. The random search (RS) conformational analysis technique using the Tripos force field was found to be capable of locating the minima found by the molecular orbital methods using systematic grid search. This suggests that the RS/Tripos force field/vacuum phase protocol is a reasonable choice for locating the local minima of MP. However, the Tripos force field gave significantly larger phenyl ring rotational barriers than the molecular orbital methods for MP and RMP. For both the neutral and protonated cases, all three methods found the phenyl ring rotational barriers for the RMP conformers/invertamers (denoted as cte, tte, and cta) to be: cte, tte > MP > cta. Solvation has negligible effect on the phenyl ring rotational barrier of RMP. The B3LYP/6-31G* density functional method was used to calculate the phenyl ring rotational barrier for neutral MP and gave results very similar to those of the HF/6-31G* method.

The mutagenic activity of 23 triazenes and, in a different set, of 24 halogenated hydroxyfuranones (MX derivatives) is quantitatively related to new features of contemporary molecular wave functions. Nowadays affordable computers are powerful enough to rapidly generate geometry-optimised ab initio wave functions at HF/3-21G*, HF/6-31G* and B3LYP/6-311 + G(2d,p) level for all molecules. The bonds of a common molecular skeleton are described by their ab initio bond lengths and local properties provided by the theory of quantum chemical topology (QCT). The chemometric analysis involves two types: one to generate a statistically validated quantitative model, and one to isolate the active center. In the former a genetic algorithm (GA) selects bond descriptors in order to optimise the cross-validation error, q2, followed by a full partial least squares (PLS) analysis, which also yields randomisation statistics. In the latter type principal components (PCs) are constructed from the original bond descriptors and their variables important to the projection (VIPs) are plotted in a histogram. This analysis suggests a preferred mechanistic pathway for the initial hydroxylation of the triazenes, an issue that has remained ambiguous so far. In the case of the hydroxyfuranones the proposed method aids the elucidation of a mechanistic ambivalence.