Benchmarking ensemble docking methods in D3R Grand Challenge 4.

The discovery of new drugs is a time consuming and expensive process. Methods such as virtual screening, which can filter out ineffective compounds from drug libraries prior to expensive experimental study, have become popular research topics. As the computational drug discovery community has grown, in order to benchmark the various advances in methodology, organizations such as the Drug Design Data Resource have begun hosting blinded grand challenges seeking to identify the best methods for ligand pose-prediction, ligand affinity ranking, and free energy calculations. Such open challenges offer a unique opportunity for researchers to partner with junior students (e.g., high school and undergraduate) to validate basic yet fundamental hypotheses considered to be uninteresting to domain experts. Here, we, a group of high school-aged students and their mentors, present the results of our participation in Grand Challenge 4 where we predicted ligand affinity rankings for the Cathepsin S protease, an important protein target for autoimmune diseases. To investigate the effect of incorporating receptor dynamics on ligand affinity rankings, we employed the Relaxed Complex Scheme, a molecular docking method paired with molecular dynamics-generated receptor conformations. We found that Cathepsin S is a difficult target for molecular docking and we explore some advanced methods such as distance-restrained docking to try to improve the correlation with experiments. This project has exemplified the capabilities of high school students when supported with a rigorous curriculum, and demonstrates the value of community-driven competitions for beginners in computational drug discovery.

Predicting Ligand Binding Kinetics Using a Markovian Milestoning with Voronoi Tessellations Multiscale Approach.

Accurate and efficient computational predictions of ligand binding kinetics can be useful to inform drug discovery campaigns, particularly in the screening and lead optimization phases. Simulation enabled estimation of kinetic rates (SEEKR) is a multiscale molecular dynamics, Brownian dynamics, and milestoning simulation approach for calculating receptor-ligand association and dissociation rates. Here, we present the implementation of a Markovian milestoning with Voronoi tessellations approach that significantly reduces the simulation cost of calculations as well as further improving their parallelizability. The new approach is applied to a host-guest system to assess its effectiveness for rank-ordering compounds by kinetic rates and to the model protein system, trypsin, with the noncovalent inhibitor benzamidine. For both applications, we demonstrate that the new approach requires up to a factor of 10 less simulation time to achieve results with comparable or increased accuracy.

Ranking of Ligand Binding Kinetics Using a Weighted Ensemble Approach and Comparison with a Multiscale Milestoning Approach.

To improve lead optimization efforts in finding the right ligand, pharmaceutical industries need to know the ligand's binding kinetics, such as binding and unbinding rate constants, which often correlate with the ligand's efficacy in vivo. To predict binding kinetics efficiently, enhanced sampling methods, such as milestoning and the weighted ensemble (WE) method, have been used in molecular dynamics (MD) simulations of these systems. However, a comparison of these enhanced sampling methods in ranking ligands has not been done. Hence, a WE approach called the concurrent adaptive sampling (CAS) algorithm that uses MD simulations was used to rank seven ligands for ß-cyclodextrin, a system in which a multiscale milestoning approach called simulation enabled estimation of kinetic rates (SEEKR) was also used, which uses both MD and Brownian dynamics simulations. Overall, the CAS algorithm can successfully rank ligands using the unbinding rate constant koff values and binding free energy ΔG values, as SEEKR did, with reduced computational cost that is about the same as SEEKR. We compare the CAS algorithm simulations with different parameters and discuss the impact of parameters in ranking ligands and obtaining rate constant and binding free energy estimates. We also discuss similarities and differences and advantages and disadvantages of SEEKR and the CAS algorithm for future use.

Multiscale simulation approaches to modeling drug-protein binding.

Simulations can provide detailed insight into the molecular processes involved in drug action, such as protein-ligand binding, and can therefore be a valuable tool for drug design and development. Processes with a large range of length and timescales may be involved, and understanding these different scales typically requires different types of simulation methodology. Ideally, simulations should be able to connect across scales, to analyze and predict how changes at one scale can influence another. Multiscale simulation methods, which combine different levels of treatment, are an emerging frontier with great potential in this area. Here we review multiscale frameworks of various types, and selected applications to biomolecular systems with a focus on drug-ligand binding.

Quantitative Ranking of Ligand Binding Kinetics with a Multiscale Milestoning Simulation Approach.

Efficient prediction and ranking of small molecule binders by their kinetic ( kon and koff) and thermodynamic ( Δ G) properties can be a valuable metric for drug lead optimization, as these quantities are often indicators of in vivo efficacy. We have previously described a hybrid molecular dynamics, Brownian dynamics, and milestoning model, Simulation Enabled Estimation of Kinetic Rates (SEEKR), that can predict kon's, koff's, and Δ G's. Here we demonstrate the effectiveness of this approach for ranking a series of seven small molecule compounds for the model system, ß-cyclodextrin, based on predicted kon's and koff's. We compare our results using SEEKR to experimentally determined rates as well as rates calculated using long time scale molecular dynamics simulations and show that SEEKR can effectively rank the compounds by koff and Δ G with reduced computational cost. We also provide a discussion of convergence properties and sensitivities of calculations with SEEKR to establish "best practices" for its future use.

SEEKR: Simulation Enabled Estimation of Kinetic Rates, A Computational Tool to Estimate Molecular Kinetics and Its Application to Trypsin-Benzamidine Binding.

We present the Simulation Enabled Estimation of Kinetic Rates (SEEKR) package, a suite of open-source scripts and tools designed to enable researchers to perform multiscale computation of the kinetics of molecular binding, unbinding, and transport using a combination of molecular dynamics, Brownian dynamics, and milestoning theory. To demonstrate its utility, we compute the kon, koff, and ΔGbind for the protein trypsin with its noncovalent binder, benzamidine, and examine the kinetics and other results generated in the context of the new software, and compare our findings to previous studies performed on the same system. We compute a kon estimate of (2.1 ± 0.3) × 107 M-1 s-1, a koff estimate of 83 ± 14 s-1, and a ΔGbind of -7.4 ± 0.1 kcal·mol-1, all of which compare closely to the experimentally measured values of 2.9 × 107 M-1 s-1, 600 ± 300 s-1, and -6.71 ± 0.05 kcal·mol-1, respectively.

Distinguishing the Protonation State of the Histidine Ligand to the Oxidized Iron-Sulfur Cluster from the MitoNEET Family of Proteins.

The iron-sulfur cluster located in the recently discovered human mitoNEET protein (and related proteins) is structurally similar to the more well-known ferredoxin and Rieske clusters. Although its biological function is uncertain, the iron-sulfur cluster in mitoNEET has been proposed to undergo proton-coupled electron transfer involving the histidine ligand to the cluster. The cluster is also released from the protein at low pH. This contribution reports density functional calculations to model the structures, vibrations, and Heisenberg coupling constants (J) for high-spin (HS), broken symmetry (BS) singlet, and extended broken symmetry (EBS) singlet states of the oxidized iron-sulfur cluster from mitoNEET. This work suggests that J values or 15 N isotopic frequency shifts may provide methods for determining experimentally whether the histidine ligand to the oxidized iron-sulfur cluster in human mitoNEET and mitoNEET-related proteins is protonated or deprotonated.

Distinguishing Protonation States of Histidine Ligands to the Oxidized Rieske Iron-Sulfur Cluster through (15) N Vibrational Frequency Shifts.

The Rieske [2Fe-2S] cluster is a vital component of many oxidoreductases, including mitochondrial cytochrome bc1; its chloroplast equivalent, cytochrome b6f; one class of dioxygenases; and arsenite oxidase. The Rieske cluster acts as an electron shuttle and its reduction is believed to couple with protonation of one of the cluster's His ligands. In cytochromes bc1 and b6f, for example, the Rieske cluster acts as the first electron acceptor in a modified Q cycle. The protonation states of the cluster's His ligands determine its ability to accept a proton and possibly an electron through a hydrogen bond to the electron carrier, ubiquinol. Experimental determination of the protonation states of a Rieske cluster's two His ligands by NMR spectroscopy is difficult, due to the close proximity of the two paramagnetic iron atoms of the cluster. Therefore, this work reports density functional calculations and proposes that difference vibrational spectroscopy with (15) N isotopic substitution may be used to assign the protonation states of the His ligands of the oxidized Rieske [2Fe-2S] complex.