Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Can We Quickly Learn to "Translate" Bioactive Molecules with Transformer Models?

Meaningful exploration of the chemical space of druglike molecules in drug design is a highly challenging task due to a combinatorial explosion of possible modifications of molecules. In this work, we address this problem with transformer models, a type of machine learning (ML) model originally developed for machine translation. By training transformer models on pairs of similar bioactive molecules from the public ChEMBL data set, we enable them to learn medicinal-chemistry-meaningful, context-dependent transformations of molecules, including those absent from the training set. By retrospective analysis on the performance of transformer models on ChEMBL subsets of ligands binding to COX2, DRD2, or HERG protein targets, we demonstrate that the models can generate structures identical or highly similar to most active ligands, despite the models having not seen any ligands active against the corresponding protein target during training. Our work demonstrates that human experts working on hit expansion in drug design can easily and quickly employ transformer models, originally developed to translate texts from one natural language to another, to "translate" from known molecules active against a given protein target to novel molecules active against the same target.

Low-data interpretable deep learning prediction of antibody viscosity using a biophysically meaningful representation.

Deep learning, aided by the availability of big data sets, has led to substantial advances across many disciplines. However, many scientific problems of practical interest lack sufficiently large datasets amenable to deep learning. Prediction of antibody viscosity is one such problem where deep learning methods have not yet been explored due to the relative scarcity of relevant training data. In this work, we overcome this limitation using a biophysically meaningful representation that enables us to develop generalizable models even under limited training data. We present, PfAbNet-viscosity, a 3D convolutional neural network architecture, to predict high-concentration viscosity of therapeutic antibodies. We show that with the electrostatic potential surface of the antibody variable region as the only input to the network, the models trained on as few as couple dozen datapoints can generalize with high accuracy. Our feature attribution analysis shows that PfAbNet-viscosity has learned key biophysical drivers of viscosity. The applicability of our approach to other biological systems is discussed.

TorsionNet: A Deep Neural Network to Rapidly Predict Small-Molecule Torsional Energy Profiles with the Accuracy of Quantum Mechanics.

Fast and accurate assessment of small-molecule dihedral energetics is crucial for molecular design and optimization in medicinal chemistry. Yet, accurate prediction of torsion energy profiles remains challenging as the current molecular mechanics (MM) methods are limited by insufficient coverage of drug-like chemical space and accurate quantum mechanical (QM) methods are too expensive. To address this limitation, we introduce TorsionNet, a deep neural network (DNN) model specifically developed to predict small-molecule torsion energy profiles with QM-level accuracy. We applied active learning to identify nearly 50k fragments (with elements H, C, N, O, F, S, and Cl) that maximized the coverage of our corporate compound library and leveraged massively parallel cloud computing resources for density functional theory (DFT) torsion scans of these fragments, generating a training data set of 1.2 million DFT energies. After training TorsionNet on this data set, we obtain a model that can rapidly predict the torsion energy profile of typical drug-like fragments with DFT-level accuracy. Importantly, our method also provides an uncertainty estimate for the predicted profiles without any additional calculations. In this report, we show that TorsionNet can accurately identify the preferred dihedral geometries observed in crystal structures. Our TorsionNet-based analysis of a diverse set of protein-ligand complexes with measured binding affinity shows a strong association between high ligand strain and low potency. We also present practical applications of TorsionNet that demonstrate how consideration of DNN-based strain energy leads to substantial improvement in existing lead discovery and design workflows. TorsionNet500, a benchmark data set comprising 500 chemically diverse fragments with DFT torsion profiles (12k MM- and DFT-optimized geometries and energies), has been created and is made publicly available.

Comprehensive Assessment of Torsional Strain in Crystal Structures of Small Molecules and Protein-Ligand Complexes using ab Initio Calculations.

The energetics of rotation around single bonds (torsions) is a key determinant of the three-dimensional shape that druglike molecules adopt in solution, the solid state, and in different biological environments, which in turn defines their unique physical and pharmacological properties. Therefore, accurate characterization of torsion angle preference and energetics is essential for the success of computational drug discovery and design. Here, we analyze torsional strain in crystal structures of druglike molecules in Cambridge structure database (CSD) and bioactive ligand conformations in protein data bank (PDB), expressing the total strain energy as a sum of strain energy from constituent rotatable bonds. We utilized cloud computing to generate torsion scan profiles of a very large collection of chemically diverse neutral fragments at DFT(B3LYP)/6-31G*//6-31G** or DFT(B3LYP)/6-31+G*//6-31+G** (for sulfur-containing molecule). With the data generated from these ab initio calculations, we performed rigorous analysis of strain due to deviation of observed torsion angles relative to their ideal gas-phase geometries. Contrary to the previous studies based on molecular mechanics, we find that in the crystalline state, molecules generally adopt low-strain conformations, with median per-torsion strain energy in CSD and PDB under one-tenth and one-third of a kcal/mol, respectively. However, for a small fraction (<5%) of motifs, external effects such as steric hindrance and hydrogen bonds result in strain penalty exceeding 2.5 kcal/mol. We find that due to poor quality of PDB structures in general, bioactive structures tend to have higher torsional strain compared to small-molecule crystal conformations. However, in the absence of structural fitting artifacts in PDB structures, protein-induced strain in bioactive conformations is quantitatively similar to those due to the packing forces in small-molecule crystal structures. This analysis allows us to establish strain energy thresholds to help identify biologically relevant conformers in a given ensemble. The work presented here is the most comprehensive study to date that demonstrates the utility and feasibility of gas-phase quantum mechanics (QM) calculations to study conformational preference and energetics of drug-size molecules. Potential applications of this study in computational lead discovery and structure-based design are discussed.

Discovery and Lead Optimization of Atropisomer D1 Agonists with Reduced Desensitization.

The discovery of D1 subtype-selective agonists with drug-like properties has been an enduring challenge for the greater part of 40 years. All known D1-selective agonists are catecholamines that bring about receptor desensitization and undergo rapid metabolism, thus limiting their utility as a therapeutic for chronic illness such as schizophrenia and Parkinson's disease. Our high-throughput screening efforts on D1 yielded a single non-catecholamine hit PF-4211 (6) that was developed into a series of potent D1 receptor agonist leads with high oral bioavailability and CNS penetration. An important structural feature of this series is the locked biaryl ring system resulting in atropisomerism. Disclosed herein is a summary of our hit-to-lead efforts on this series of D1 activators culminating in the discovery of atropisomer 31 (PF-06256142), a potent and selective orthosteric agonist of the D1 receptor that has reduced receptor desensitization relative to dopamine and other catechol-containing agonists.

Fast and accurate generation of ab initio quality atomic charges using nonparametric statistical regression.

We introduce a class of partial atomic charge assignment method that provides ab initio quality description of the electrostatics of bioorganic molecules. The method uses a set of models that neither have a fixed functional form nor require a fixed set of parameters, and therefore are capable of capturing the complexities of the charge distribution in great detail. Random Forest regression is used to build separate charge models for elements H, C, N, O, F, S, and Cl, using training data consisting of partial charges along with a description of their surrounding chemical environments; training set charges are generated by fitting to the b3lyp/6-31G* electrostatic potential (ESP) and are subsequently refined to improve consistency and transferability of the charge assignments. Using a set of 210 neutral, small organic molecules, the absolute hydration free energy calculated using these charges in conjunction with Generalized Born solvation model shows a low mean unsigned error, close to 1 kcal/mol, from the experimental data. Using another large and independent test set of chemically diverse organic molecules, the method is shown to accurately reproduce charge-dependent observables--ESP and dipole moment--from ab initio calculations. The method presented here automatically provides an estimate of potential errors in the charge assignment, enabling systematic improvement of these models using additional data. This work has implications not only for the future development of charge models but also in developing methods to describe many other chemical properties that require accurate representation of the electronic structure of the system.

Exploring aromatic chemical space with NEAT: novel and electronically equivalent aromatic template.

In this paper, we describe a lead transformation tool, NEAT (Novel and Electronically equivalent Aromatic Template), which can help identify novel aromatic rings that are estimated to have similar electrostatic potentials, dipoles, and hydrogen bonding capabilities to a query template; hence, they may offer similar bioactivity profiles. In this work, we built a comprehensive heteroaryl database, and precalculated high-level quantum mechanical (QM) properties, including electrostatic potential charges, hydrogen bonding ability, dipole moments, chemical reactivity, and othe properties. NEAT bioisosteric similarities are based on the electrostatic potential surface calculated by Brood, using the precalculated QM ESP charges and other QM properties. Compared with existing commercial lead transformation software, (1) NEAT is the only one that covers the comprehensive heteroaryl chemical space, and (2) NEAT offers a better characterization of novel aryl cores by using high-evel QM properties that are relevant to molecular interactions. NEAT provides unique value to medicinal chemists quickly exploring the largely uncharted aromatic chemical space, and one successful example of its application is discussed herein.

Structure-based druggability assessment--identifying suitable targets for small molecule therapeutics.

A target is druggable if it can be modulated in vivo by a drug-like molecule. The general properties of oral drugs are summarized by the 'rule of 5' which specifies parameters related to size and lipophilicity. Structure-based target druggability assessment consists of predicting ligand-binding sites on the protein that are complementary to these drug-like properties. Automated identification of ligand-binding sites can use geometrical considerations alone or include specific physicochemical properties of the protein surface. Features of a pocket's size and shape, together with measures of its hydrophobicity, are most informative in identifying suitable drug-binding pockets. The recent availability of several validation sets of druggable versus undruggable targets has helped fuel the development of more elaborate methods.

Understanding the impact of the P-loop conformation on kinase selectivity.

This work addresses the link between selectivity and an unusual, folded conformation for the P-loop observed initially for MAP4K4 and subsequently for other kinases. Statistical and computational analyses of our crystal structure database demonstrate that inhibitors that induce the P-loop folded conformation tend to be more selective, especially if they take advantage of this specific conformation by interacting more favorably with a conserved Tyr or Phe residue from the P-loop.

Modeling G protein-coupled receptors for structure-based drug discovery using low-frequency normal modes for refinement of homology models: application to H3 antagonists.

G Protein-Coupled Receptors (GPCRs) are integral membrane proteins that play important role in regulating key physiological functions, and are targets of about 50% of all recently launched drugs. High-resolution experimental structures are available only for very few GPCRs. As a result, structure-based drug design efforts for GPCRs continue to rely on in silico modeling, which is considered to be an extremely difficult task especially for these receptors. Here, we describe Gmodel, a novel approach for building 3D atomic models of GPCRs using a normal mode-based refinement of homology models. Gmodel uses a small set of relevant low-frequency vibrational modes derived from Random Elastic Network model to efficiently sample the large-scale receptor conformation changes and generate an ensemble of alternative models. These are used to assemble receptor-ligand complexes by docking a known active into each of the alternative models. Each of these is next filtered using restraints derived from known mutation and binding affinity data and is refined in the presence of the active ligand. In this study, Gmodel was applied to generate models of the antagonist form of histamine 3 (H3) receptor. The validity of this novel modeling approach is demonstrated by performing virtual screening (using the refined models) that consistently produces highly enriched hit lists. The models are further validated by analyzing the available SAR related to classical H3 antagonists, and are found to be in good agreement with the available experimental data, thus providing novel insights into the receptor-ligand interactions.

Comparison of several molecular docking programs: pose prediction and virtual screening accuracy.

Molecular docking programs are widely used modeling tools for predicting ligand binding modes and structure based virtual screening. In this study, six molecular docking programs (DOCK, FlexX, GLIDE, ICM, PhDOCK, and Surflex) were evaluated using metrics intended to assess docking pose and virtual screening accuracy. Cognate ligand docking to 68 diverse, high-resolution X-ray complexes revealed that ICM, GLIDE, and Surflex generated ligand poses close to the X-ray conformation more often than the other docking programs. GLIDE and Surflex also outperformed the other docking programs when used for virtual screening, based on mean ROC AUC and ROC enrichment values obtained for the 40 protein targets in the Directory of Useful Decoys (DUD). Further analysis uncovered general trends in accuracy that are specific for particular protein families. Modifying basic parameters in the software was shown to have a significant effect on docking and virtual screening results, suggesting that expert knowledge is critical for optimizing the accuracy of these methods.

Comparative protein structure modeling by combining multiple templates and optimizing sequence-to-structure alignments.

MOTIVATION: Two major bottlenecks in advancing comparative protein structure modeling are the efficient combination of multiple template structures and the generation of a correct input target-template alignment. RESULTS: A novel method, Multiple Mapping Method with Multiple Templates (M4T) is introduced that implements an algorithm to automatically select and combine Multiple Template structures (MT) and an alignment optimization protocol (Multiple Mapping Method, MMM). The MT module of M4T selects and combines multiple template structures through an iterative clustering approach that takes into account the 'unique' contribution of each template, their sequence similarity among themselves and to the target sequence, and their experimental resolution. MMM is a sequence-to-structure alignment method that optimally combines alternatively aligned regions according to their fit in the structural environment of the template structure. The resulting M4T alignment is used as input to a comparative modeling module. The performance of M4T has been benchmarked on CASP6 comparative modeling target sequences and on a larger independent test set, and showed favorable performance to current state of the art methods.

MMM: a sequence-to-structure alignment protocol.

MOTIVATION: Accurate alignment of a target sequence to a template structure continues to be a bottleneck in producing good quality comparative protein structure models. RESULTS: Multiple Mapping Method (MMM) is a comparative protein structure modeling server with an emphasis on a novel alignment optimization protocol. MMM takes inputs from five profile-to-profile based alignment methods. The alternatively aligned regions from the input alignment set are combined according to their fit in the structural environment of the template structure. The resulting, optimally spliced MMM alignment is used as input to an automated comparative modeling module to produce a full atom model. AVAILABILITY: The MMM server is freely accessible at http://www.fiserlab.org/servers/mmm

Multiple mapping method: a novel approach to the sequence-to-structure alignment problem in comparative protein structure modeling.

A major bottleneck in comparative protein structure modeling is the quality of input alignment between the target sequence and the template structure. A number of alignment methods are available, but none of these techniques produce consistently good solutions for all cases. Alignments produced by alternative methods may be superior in certain segments but inferior in others when compared to each other; therefore, an accurate solution often requires an optimal combination of them. To address this problem, we have developed a new approach, Multiple Mapping Method (MMM). The algorithm first identifies the alternatively aligned regions from a set of input alignments. These alternatively aligned segments are scored using a composite scoring function, which determines their fitness within the structural environment of the template. The best scoring regions from a set of alternative segments are combined with the core part of the alignments to produce the final MMM alignment. The algorithm was tested on a dataset of 1400 protein pairs using 11 combinations of two to four alignment methods. In all cases MMM showed statistically significant improvement by reducing alignment errors in the range of 3 to 17%. MMM also compared favorably over two alignment meta-servers. The algorithm is computationally efficient; therefore, it is a suitable tool for genome scale modeling studies.

Single-atom test of all-atom empirical potentials: Fe in myoglobin.

The measured Fe vibrational density of states in deoxy-myoglobin, obtained from nuclear resonance vibrational spectroscopy, is compared to results from a normal-mode analysis using an all-atom empirical potential. Substantial disagreement reveals that for this one atom, the empirical potential does not accurately describe the actual forces. A Green function technique is developed to calculate the iron vibrational spectrum of deoxy-myoglobin by coupling the independently calculated heme and globin normal modes; nonbonded interactions between the heme molecule and the protein are essential for a good fit to the measurements. A projection of the eigenvectors from this potential onto the displacements induced by binding of CO demonstrates that normal modes over a broad range centered around 50-150 cm(-1) may drive the ligand-induced structural changes.

Direct determination of the complete set of iron normal modes in a porphyrin-imidazole model for carbonmonoxy-heme proteins: [Fe(TPP)(CO)(1-MeIm)].

Detailed Fe vibrational spectra have been obtained for the heme model complex [Fe(TPP)(CO)(1-MeIm)] using a new, highly selective and quantitative technique, Nuclear Resonance Vibrational Spectroscopy (NRVS). This spectroscopy measures the complete vibrational density of states for iron atoms, from which normal modes can be calculated via refinement of the force constants. These data and mode assignments can reveal previously undetected vibrations and are useful for validating predictions based on optical spectroscopies and density functional theory, for example. Vibrational modes of the iron porphyrin-imidazole compound [Fe(TPP)(CO)(1-MeIm)] have been determined by refining normal mode calculations to NRVS data obtained at an X-ray synchrotron source. Iron dynamics of this compound, which serves as a useful model for the active site in the six-coordinate heme protein, carbonmonoxy-myoglobin, are discussed in relation to recently determined dynamics of a five-coordinate deoxy-myoglobin model, [Fe(TPP)(2-MeHIm)]. For the first time in a six-coordinate heme system, the iron-imidazole stretch mode has been observed, at 226 cm(-)(1). The heme in-plane modes with large contributions from the nu(42), nu(49), nu(50), and nu(53) modes of the core porphyrin are identified. In general, the iron modes can be attributed to coupling with the porphyrin core, the CO ligand, the imidazole ring, and/or the phenyl rings. Other significant findings are the observation that the porphyrin ring peripheral substituents are strongly coupled to the iron doming mode and that the Fe-C-O tilting and bending modes are related by a negative interaction force constant.

Iron normal mode dynamics in (nitrosyl)iron(II)tetraphenylporphyrin from X-ray nuclear resonance data.

The complete iron atom vibrational spectrum has been obtained by refinement of normal mode calculations to nuclear inelastic x-ray absorption data from (nitrosyl)iron(II)tetraphenylporphyrin, FeTPP(NO), a useful model for heme dynamics in myoglobin and other heme proteins. Nuclear resonance vibrational spectroscopy (NRVS) provides a direct measurement of the frequency and iron amplitude for all normal modes involving significant displacement of (57)Fe. The NRVS measurements on isotopically enriched single crystals permit determination of heme in-plane and out-of-plane modes. Excellent agreement between the calculated and experimental values of frequency and iron amplitude for each mode is achieved by a force-field refinement. Significantly, we find that the presence of the phenyl groups and the NO ligand leads to substantial mixing of the porphyrin core modes. This first picture of the entire iron vibrational density of states for a porphyrin compound provides an improved model for the role of iron atom dynamics in the biological functioning of heme proteins.

Iron normal mode dynamics in a porphyrin-imidazole model for deoxyheme proteins.

Iron vibrational modes of a deoxyheme protein model (2-methylimidazole)(tetraphenylporphinato)iron(II), [Fe(TPP)(2-MeImH)], have been studied by refining normal mode calculations to nuclear resonance vibrational spectroscopy (NRVS) data. The NRVS measurements give quantitative frequencies and iron amplitudes of all modes with significant Fe vibrational motion. Modes with in-plane displacement of iron are distinguished from those involving out-of-plane motion by measurements on oriented single-crystal samples. Normal modes having large overlaps with in-plane nu(42), nu(50), and nu(53) modes of the porphyrin core are identified, as well as several modes with large iron-imidazole stretch components. An out-of-plane mode at 78 cm(-1) shows significant doming of the porphyrin core, but the largest Fe doming motion arises from the coupling of phenyls and imidazole at 25 cm(-1).