AutoDesigner - Core Design, a De Novo Design Algorithm for Chemical Scaffolds: Application to the Design and Synthesis of Novel Selective Wee1 Inhibitors.

The hit identification stage of a drug discovery program generally involves the design of novel chemical scaffolds with desired biological activity against the target(s) of interest. One common approach is scaffold hopping, which is the manual design of novel scaffolds based on known chemical matter. One major limitation of this approach is narrow chemical space exploration, which can lead to difficulties in maintaining or improving biological activity, selectivity, and favorable property space. Another limitation is the lack of preliminary structure-activity relationship (SAR) data around these designs, which could lead to selecting suboptimal scaffolds to advance lead optimization. To address these limitations, we propose AutoDesigner - Core Design (CoreDesign), a de novo scaffold design algorithm. Our approach is a cloud-integrated, de novo design algorithm for systematically exploring and refining chemical scaffolds against biological targets of interest. The algorithm designs, evaluates, and optimizes a vast range, from millions to billions, of molecules in silico, following defined project parameters encompassing structural novelty, physicochemical attributes, potency, and selectivity using active-learning FEP. To validate CoreDesign in a real-world drug discovery setting, we applied it to the design of novel, potent Wee1 inhibitors with improved selectivity over PLK1. Starting from a single known ligand and receptor structure, CoreDesign rapidly explored over 23 billion molecules to identify 1,342 novel chemical series with a mean of 4 compounds per scaffold. To rapidly analyze this large amount of data and prioritize chemical scaffolds for synthesis, we utilize t-Distributed Stochastic Neighbor Embedding (t-SNE) plots of in silico properties. The chemical space projections allowed us to rapidly identify a structurally novel 5-5 fused core meeting all the hit-identification requirements. Several compounds were synthesized and assayed from the scaffold, displaying good potency against Wee1 and excellent PLK1 selectivity. Our results suggest that CoreDesign can significantly speed up the hit-identification process and increase the probability of success of drug discovery campaigns by allowing teams to bring forward high-quality chemical scaffolds derisked by the availability of preliminary SAR.

Exploiting high-energy hydration sites for the discovery of potent peptide aldehyde inhibitors of the SARS-CoV-2 main protease with cellular antiviral activity.

Small-molecule antivirals that prevent the replication of the SARS-CoV-2 virus by blocking the enzymatic activity of its main protease (Mpro) are and will be a tenet of pandemic preparedness. However, the peptidic nature of such compounds often precludes the design of compounds within favorable physical property ranges, limiting cellular activity. Here we describe the discovery of peptide aldehyde Mpro inhibitors with potent enzymatic and cellular antiviral activity. This structure-activity relationship (SAR) exploration was guided by the use of calculated hydration site thermodynamic maps (WaterMap) to drive potency via displacement of waters from high-energy sites. Thousands of diverse compounds were designed to target these high-energy hydration sites and then prioritized for synthesis by physics- and structure-based Free-Energy Perturbation (FEP+) simulations, which accurately predicted biochemical potencies. This approach ultimately led to the rapid discovery of lead compounds with unique SAR that exhibited potent enzymatic and cellular activity with excellent pan-coronavirus coverage.

OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins.

The parametrization and validation of the OPLS3 force field for small molecules and proteins are reported. Enhancements with respect to the previous version (OPLS2.1) include the addition of off-atom charge sites to represent halogen bonding and aryl nitrogen lone pairs as well as a complete refit of peptide dihedral parameters to better model the native structure of proteins. To adequately cover medicinal chemical space, OPLS3 employs over an order of magnitude more reference data and associated parameter types relative to other commonly used small molecule force fields (e.g., MMFF and OPLS_2005). As a consequence, OPLS3 achieves a high level of accuracy across performance benchmarks that assess small molecule conformational propensities and solvation. The newly fitted peptide dihedrals lead to significant improvements in the representation of secondary structure elements in simulated peptides and native structure stability over a number of proteins. Together, the improvements made to both the small molecule and protein force field lead to a high level of accuracy in predicting protein-ligand binding measured over a wide range of targets and ligands (less than 1 kcal/mol RMS error) representing a 30% improvement over earlier variants of the OPLS force field.

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field.

Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

Leveraging Data Fusion Strategies in Multireceptor Lead Optimization MM/GBSA End-Point Methods.

Accurate and efficient affinity calculations are critical to enhancing the contribution of in silico modeling during the lead optimization phase of a drug discovery campaign. Here, we present a large-scale study of the efficacy of data fusion strategies to leverage results from end-point MM/GBSA calculations in multiple receptors to identify potent inhibitors among an ensemble of congeneric ligands. The retrospective analysis of 13 congeneric ligand series curated from publicly available data across seven biological targets demonstrates that in 90% of the individual receptor structures MM/GBSA scores successfully identify subsets of inhibitors that are more potent than a random selection, and data fusion strategies that combine MM/GBSA scores from each of the receptors significantly increase the robustness of the predictions. Among nine different data fusion metrics based on consensus scores or receptor rankings, the SumZScore (i.e., converting MM/GBSA scores into standardized Z-Scores within a receptor and computing the sum of the Z-Scores for a given ligand across the ensemble of receptors) is found to be a robust and physically meaningful metric for combining results across multiple receptors. Perhaps most surprisingly, even with relatively low to modest overall correlations between SumZScore and experimental binding affinities, SumZScore tends to reliably prioritize subsets of inhibitors that are at least as potent as those that are prioritized from a "best" single receptor identified from known compounds within the congeneric series.

Deconstructing activation events in rhodopsin.

Activation of class-A G-protein-coupled receptors (GPCRs) involves large-scale reorganization of the H3/H6 interhelical network. In rhodopsin (Rh), this process is coupled to a change in the protonation state of a key residue, E134, whose exact role in activation is not well understood. Capturing this millisecond pH-dependent process is a well-appreciated challenge. We have developed a scheme combining the harmonic Fourier beads (HFB) method and constant-pH molecular dynamics with pH-based replica exchange (pH-REX) to gain insight into the structural changes that occur along the activation pathway as a function of the protonation state of E134. Our results indicate that E134 is protonated as a consequence of tilting of H6 by ca. 4.0° with respect to its initial position and simultaneous rotation by ca. 23° along its principal axis. The movement of H6 is associated with breakage of the E247-R135 and R135-E134 salt bridges and concomitant release of the E134 side chain, which results in an increase in its pKa value above physiological pH. An increase in the hydrophobicity of the environment surrounding E134 leads to further tilting and rotation of H6 and upshift of the E134 pKa. Such atomic-level information, which is not accessible through experiments, refines the earlier proposed sequential model of Rh activation (see: Zaitseva, E.; et al. Sequential Rearrangement of Interhelical Networks Upon Rhodopsin Activation in Membranes: The Meta IIa Conformational Substate . J. Am. Chem. Soc. 2010, 132, 4815) and argues that the E134 protonation switch is both a cause and a consequence of the H6 motion.

Towards Accurate Prediction of Protonation Equilibrium of Nucleic Acids.

The role of protonated nucleotides in modulating the pH-dependent properties of nucleic acids is one of the emerging frontiers in the field of nucleic acid biology. The recent development of a constant pH molecular dynamics simulation (CPHMDMSλD) framework for simulating nucleic acids has provided a tool for realistic simulations of pH-dependent dynamics. We enhanced the CPHMDMSλD framework with pH-based replica exchange (pH-REX), which significantly improves the sampling of both titration and spatial coordinates. The results from our pKa calculations for the GAAA tetraloop, which was predicted with lower accuracy previously due to sampling challenges, demonstrates that pH-REX reduces the average unsigned error (AUE) to 0.7 pKa units, and the error of the most poorly predicted residue A17 was drastically reduced from 2.9 to 1.2 pKa unit. Lastly, we show that pH-REX CPHMDMSλD simulations can be used to identify the dominant conformation of nucleic acid structures in alternate pH environments. This work suggests that pH-REX CPHMDMSλD simulations provide a practical tool for predicting nucleic acid protonation equilibrium from first-principles, and offering structural and mechanistic insight into the study of pH-dependent properties of nucleic acids.

pH-dependent dynamics of complex RNA macromolecules.

The role of pH-dependent protonation equilibrium in modulating RNA dynamics and function is one of the key unanswered questions in RNA biology. Molecular dynamics (MD) simulations can provide insight into the mechanistic roles of protonated nucleotides, but it is only capable of modeling fixed protonation states and requires prior knowledge of the key residue's protonation state. Recently, we developed a framework for constant pH molecular dynamics simulations (CPHMDMSλD) of nucleic acids, where the nucleotides' protonation states are modeled as dynamic variables that are coupled to the structural dynamics of the RNA. In the present study, we demonstrate the application of CPHMDMSλD to the lead-dependent ribozyme; establishing the validity of this approach for modeling complex RNA structures. We show that CPHMDMSλD accurately predicts the direction of the pKa shifts and reproduces experimentally-measured microscopic pKa values with an average unsigned error of 1.3 pKa units. The effects of coupled titration states in RNA structures are modeled, and the importance of conformation sampling is highlighted. The general accuracy of CPHMDMSλD simulations in reproducing pH-dependent observables reported in this work demonstrates that constant pH simulations provides a powerful tool to investigate pH-dependent processes in nucleic acids.

Assessing the quality of absolute hydration free energies among CHARMM-compatible ligand parameterization schemes.

Multipurpose atom-typer for CHARMM (MATCH), an atom-typing toolset for molecular mechanics force fields, was recently developed in our laboratory. Here, we assess the ability of MATCH-generated parameters and partial atomic charges to reproduce experimental absolute hydration free energies for a series of 457 small neutral molecules in GBMV2, Generalized Born with a smooth SWitching (GBSW), and fast analytical continuum treatment of solvation (FACTS) implicit solvent models. The quality of hydration free energies associated with small molecule parameters obtained from ParamChem, SwissParam, and Antechamber are compared. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, these automated parameterization schemes with GBMV2 and GBSW demonstrate reasonable agreement with experimental hydration free energies (average unsigned errors of 0.9-1.5 kcal/mol and R(2) of 0.63-0.87). GBMV2 and GBSW consistently provide slightly more accurate estimates than FACTS, whereas Antechamber parameters yield marginally more accurate estimates than the current generation of MATCH, ParamChem, and SwissParam parameterization strategies. Modeling with MATCH libraries that are derived from different CHARMM topology and parameter files highlights the importance of having sufficient coverage of chemical space within the underlying databases of these automated schemes and the benefit of targeting specific functional groups for parameterization efforts to maximize both the breadth and the depth of the parameterized space.

Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors.

Accurate and reliable calculation of protein-ligand binding affinities remains a hotbed of computer-aided drug design research. Despite the potentially large impact FEP (free energy perturbation) may have in drug design projects, practical applications of FEP in industrial contexts have been limited. In this work, we use a recently developed method, FEP/REST (free energy perturbation/replica exchange with solute tempering), to calculate the relative binding affinities for a set of congeneric ligands binding to the CDK2 receptor. We compare the FEP/REST results with traditional FEP/MD (molecular dynamics) results and MM/GBSA (molecular mechanics/Generalized Born Surface Area model) results and examine why FEP/REST performed notably better than these other methods, as well as why certain ligand mutations lead to large increases of the binding affinity while others do not. We also introduce a mathematical framework for assessing the consistency and reliability of the calculations using cycle closures in FEP mutation paths.

Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent.

The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, respectively, and are assumed to be unprotonated under physiological conditions. However, evidence from recent NMR and X-Ray crystallography studies has revealed the prevalence of protonated adenine and cytosine in RNA macromolecules. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first constant pH MD simulations (CPHMD) for nucleic acids in explicit solvent in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λ(Nexp) for λ that was recently developed for Multi-Site λ-Dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calculated pKa values of simple nucleotides are in a good agreement with experimentally measured values, with a mean absolute error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biological properties of RNA macromolecules.

MATCH: an atom-typing toolset for molecular mechanics force fields.

We introduce a toolset of program libraries collectively titled multipurpose atom-typer for CHARMM (MATCH) for the automated assignment of atom types and force field parameters for molecular mechanics simulation of organic molecules. The toolset includes utilities for the conversion of multiple chemical structure file formats into a molecular graph. A general chemical pattern-matching engine using this graph has been implemented whereby assignment of molecular mechanics atom types, charges, and force field parameters are achieved by comparison against a customizable list of chemical fragments. While initially designed to complement the CHARMM simulation package and force fields by generating the necessary input topology and atom-type data files, MATCH can be expanded to any force field and program, and has core functionality that makes it extendable to other applications such as fragment-based property prediction. In this work, we demonstrate the accurate construction of atomic parameters of molecules within each force field included in CHARMM36 through exhaustive cross validation studies illustrating that bond charge increment rules derived from one force field can be transferred to another. In addition, using leave-one-out substitution it is shown that it is also possible to substitute missing intra and intermolecular parameters with ones included in a force field to complete the parameterization of novel molecules. Finally, to demonstrate the robustness of MATCH and the coverage of chemical space offered by the recent CHARMM general force field (Vanommeslaeghe, et al., J Comput Chem 2010, 31, 671), one million molecules from the PubChem database of small molecules are typed, parameterized, and minimized.

Multi-Site λ-dynamics for simulated Structure-Activity Relationship studies.

Multi-Site λ-dynamics (MSλD) is a new free energy simulation method that is based on λ-dynamics. It has been developed to enable multiple substituents at multiple sites on a common ligand core to be modeled simultaneously and their free energies assessed. The efficacy of MSλD for estimating relative hydration free energies and relative binding affinties is demonstrated using three test systems. Model compounds representing multiple identical benzene, dihydroxybenzene and dimethoxybenzene molecules show total combined MSλD trajectory lengths of ~1.5 ns are sufficient to reliably achieve relative hydration free energy estimates within 0.2 kcal/mol and are less sensitive to the number of trajectories that are used to generate these estimates for hybrid ligands that contain up to ten substituents modeled at a single site or five substituents modeled at each of two sites. Relative hydration free energies among six benzene derivatives calculated from MSλD simulations are in very good agreement with those from alchemical free energy simulations (with average unsigned differences of 0.23 kcal/mol and R(2)=0.991) and experiment (with average unsigned errors of 1.8 kcal/mol and R(2)=0.959). Estimates of the relative binding affinities among 14 inhibitors of HIV-1 reverse transcriptase obtained from MSλD simulations are in reasonable agreement with those from traditional free energy simulations and experiment (average unsigned errors of 0.9 kcal/mol and R(2)=0.402). For the same level of accuracy and precision MSλD simulations are achieved ~20-50 times faster than traditional free energy simulations and thus with reliable force field parameters can be used effectively to screen tens to hundreds of compounds in structure-based drug design applications.

Applying efficient implicit nongeometric constraints in alchemical free energy simulations.

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.

Surveying implicit solvent models for estimating small molecule absolute hydration free energies.

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1-BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0-1.7 kcal/mol and R(2)=0.81-0.91) and with experimental hydration free energies (average unsigned errors=1.1-1.4 kcal/mol and R(2)=0.66-0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.

On the mechanism of crystalline polymorph selection by polymer heteronuclei.

The phase-selective crystallization of acetaminophen (ACM) using insoluble polymers as heteronuclei was investigated in a combined experimental and computational effort to elucidate the mechanism of polymer-induced heteronucleation (PIHn). ACM heteronucleates from supersaturated aqueous solution in its most thermodynamically stable monoclinic form on poly(n-butyl methacrylate), whereas the metastable orthorhombic form is observed on poly(methyl methacrylate). When ACM crystals were grown through vapor deposition, only the monoclinic polymorph was observed on each polymer. Each crystallization condition leads to a unique powder X-ray diffraction pattern with the major preferred orientation corresponding to the crystallographic faces in which these crystal phases nucleate from surfaces of the polymers. The molecular recognition events leading to these outcomes are elucidated with the aid of computed polymer-crystal binding energies using docking simulations. This investigation illuminates the mechanism by which phase selection occurs during the crystallization of ACM using polymers as heteronuclei, paving the way for the improvement of methods for polymorph selection and discovery based on heterogeneous nucleation promoters.

Lambda-dynamics free energy simulation methods.

Free energy calculations are fundamental to obtaining accurate theoretical estimates of many important biological phenomena including hydration energies, protein-ligand binding affinities and energetics of conformational changes. Unlike traditional free energy perturbation and thermodynamic integration methods, lambda-dynamics treats the conventional "lambda" as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present article, we provide an overview of the theory of lambda-dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how lambda-dynamics has been used to rapidly and reliably compute relative hydration free energies and binding affinities for series of ligands, to accurately identify crystallographically observed binding modes starting from incorrect orientations, and to model the effects of mutations upon protein stability. Finally, we suggest how lambda-dynamics may be extended to facilitate modeling efforts in structure-based drug design.

Validating CHARMM parameters and exploring charge distribution rules in structure-based drug design.

Using an extensive series of TIBO compounds that are non-nucleoside inhibitors of HIV-1 reverse transcriptase, we have systematically evaluated the quality of recently developed ligand parameters that are consistent with the CHARMM22 force field. Thermodynamic integration simulations for 44 pairs of TIBO compounds achieve a high level of success with an overall average unsigned error (AUE) in the relative binding affinities of 1.3 kcal/mol; however, the accuracy is strongly dependent on the size differential between the substituents sampled as well as the class of functional group. Low errors are observed among the alkyl, allyl, aldehyde, nitrile, trifluorinated methyl, and halide TIBO derivatives and large systematic errors among thioether derivatives. We have also investigated how different charge assignment schemes for small molecules impact the quality of computed binding affinities for a subset of this series. This study demonstrates the advantage of using model compounds to derive physically meaningful charge distributions and bond-charge increments for rapidly expanding fragment libraries for drug development applications. Specifically, in the absence of a bond-charge increment for a given pair of atom types, the strategy of adopting CHELPG charges from localized regions of model compounds provides reliable results when modeling with the CHARMM force field.

Exploring structural variability in X-ray crystallographic models using protein local optimization by torsion-angle sampling.

Modeling structural variability is critical for understanding protein function and for modeling reliable targets for in silico docking experiments. Because of the time-intensive nature of manual X-ray crystallographic refinement, automated refinement methods that thoroughly explore conformational space are essential for the systematic construction of structurally variable models. Using five proteins spanning resolutions of 1.0-2.8 A, it is demonstrated how torsion-angle sampling of backbone and side-chain libraries with filtering against both the chemical energy, using a modern effective potential, and the electron density, coupled with minimization of a reciprocal-space X-ray target function, can generate multiple structurally variable models which fit the X-ray data well. Torsion-angle sampling as implemented in the Protein Local Optimization Program (PLOP) has been used in this work. Models with the lowest R(free) values are obtained when electrostatic and implicit solvation terms are included in the effective potential. HIV-1 protease, calmodulin and SUMO-conjugating enzyme illustrate how variability in the ensemble of structures captures structural variability that is observed across multiple crystal structures and is linked to functional flexibility at hinge regions and binding interfaces. An ensemble-refinement procedure is proposed to differentiate between variability that is a consequence of physical conformational heterogeneity and that which reflects uncertainty in the atomic coordinates.

HIV-1 reverse transcriptase structure with RNase H inhibitor dihydroxy benzoyl naphthyl hydrazone bound at a novel site.

The rapid emergence of drug-resistant variants of human immunodeficiency virus, type 1 (HIV-1), has limited the efficacy of anti-acquired immune deficiency syndrome (AIDS) treatments, and new lead compounds that target novel binding sites are needed. We have determined the 3.15 A resolution crystal structure of HIV-1 reverse transcriptase (RT) complexed with dihydroxy benzoyl naphthyl hydrazone (DHBNH), an HIV-1 RT RNase H (RNH) inhibitor (RNHI). DHBNH is effective against a variety of drug-resistant HIV-1 RT mutants. While DHBNH has little effect on most aspects of RT-catalyzed DNA synthesis, at relatively high concentrations it does inhibit the initiation of RNA-primed DNA synthesis. Although primarily an RNHI, DHBNH binds >50 A away from the RNH active site, at a novel site near both the polymerase active site and the non-nucleoside RT inhibitor (NNRTI) binding pocket. When DHBNH binds, both Tyr181 and Tyr188 remain in the conformations seen in unliganded HIV-1 RT. DHBNH interacts with conserved residues (Asp186, Trp229) and has substantial interactions with the backbones of several less well-conserved residues. On the basis of this structure, we designed substituted DHBNH derivatives that interact with the NNRTI-binding pocket. These compounds inhibit both the polymerase and RNH activities of RT.