Milestoning simulation of ligand dissociation from the glycogen synthase kinase 3ß.

As drug-binding kinetics has become an important factor to be considered in modern drug discovery, this work evaluated the ability of the Milestoning method in computing the absolute dissociation rate of a ligand from the serine-threonine kinase, glycogen synthase kinase 3ß, which is a target for designing drugs to treat diseases such as neurodegenerative disorders and diabetes. We found that the Milestoning method gave good agreement with experiment with modest computational costs. Although the time scale for dissociation lasted tens of seconds, the collective molecular dynamics simulations total less than 1µs. Computing the committor function helped to identify the transition states (TSs), in which the ligand moved substantially away from the binding pocket. The glycine-rich loop with a serine residue attaching to its tips was found to undergo large movement from the bound to the TSs and might play a role in controlling drug-dissociation kinetics.

How are N-methylcarbamates encapsulated by ß-cyclodextrin: insight into the binding mechanism.

Guest molecules containing chromophore groups encapsulated by ß-cyclodextrin (ß-CD) generate circular dichroism (CD) signals, which enables a preliminary prediction of their binding modes. However, the accurate determination of the representative binding conformation (RC) remains a challenging task due to the complex conformational space of these host-guest systems. Here, we combine a molecular dynamics/quantum mechanics/continuum solvent model (MD/QM/CSM) with induced circular dichroism (ICD) data (N. L. Pacioni, A. B. Pierini and A. V. Veglia, Spectrochim. Acta A Mol. Biomol. Spectrosc., 2013, 103, 319-324.) to explore the binding mechanism of ß-CD with four N-methylcarbamate molecules: promecarb (PC), bendiocarb (BC), carbaryl (CY) and carbofuran (CF). In aqueous solution, their stability decreases as: PC > BC > CY > CF. Comparing the ECD spectra computed from TD-DFT with the ICD data can help eliminate many common binding configurations and identify the RC. The host-guest binding affinities (BAs) estimated using a ONIOM2(B971:PM6)/SMD model reproduce the measured binding trend, reveal the competition between the non-covalent interaction and solvent effect and explain the large difference in their binding modes. We also examine the fluctuations in the computed BA using similar structures.

Simulation of ligand dissociation kinetics from the protein kinase PYK2.

Early-stage drug discovery projects often focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. The kinetics of drug binding are ignored but can have significant influence on drug efficacy. Therefore, increasing attention has been paid on evaluating drug-binding kinetics early in a drug discovery process. Simulating drug-binding kinetics at the atomic level is challenging for the long time scale involved. Here, we used the transition-based reweighting analysis method (TRAM) with the Markov state model to study the dissociation of a ligand from the protein kinase PYK2. TRAM combines biased and unbiased simulations to reduce computational costs. This work used the umbrella sampling technique for the biased simulations. Although using the potential of mean force from umbrella sampling simulations with the transition-state theory over-estimated the dissociation rate by three orders of magnitude, TRAM gave a dissociation rate within an order of magnitude of the experimental value.

Using machine learning to improve ensemble docking for drug discovery.

Ensemble docking has provided an inexpensive method to account for receptor flexibility in molecular docking for virtual screening. Unfortunately, as there is no rigorous theory to connect the docking scores from multiple structures to measured activity, researchers have not yet come up with effective ways to use these scores to classify compounds into actives and inactives. This shortcoming has led to the decrease, rather than an increase in the performance of classifying compounds when more structures are added to the ensemble. Previously, we suggested machine learning, implemented in the form of a naïve Bayesian model could alleviate this problem. However, the naïve Bayesian model assumed that the probabilities of observing the docking scores to different structures to be independent. This approximation might prevent it from achieving even higher performance. In the work presented in this paper, we have relaxed this approximation when using several other machine learning methods-k nearest neighbor, logistic regression, support vector machine, and random forest-to improve ensemble docking. We found significant improvement.

Steered molecular dynamics simulations for uncovering the molecular mechanisms of drug dissociation and for drug screening: A test on the focal adhesion kinase.

Drug-binding kinetics could play important roles in determining the efficacy of drugs and has caught the attention of more drug designers. Using the dissociation of 1H-pyrrolo[2,3-b]-pyridines from the focal adhesion kinase as an example, this work finds that steered molecular dynamics simulations could help screen compounds with long-residence times. It also reveals a two-step mechanism of ligand dissociation resembling the release of ADP from protein kinase A reported earlier. A phenyl group attaching to the pyrrole prolongs residence time by creating a large activation barrier for transition from the bound to the intermediate state when it becomes exposed to the solvent. Principal component analysis shows that ligand dissociation does not couple with large-scale collective motions of the protein involving many of its amino acids. Rather, a small subset of amino acids dominates. Some of these amino acids do not contact the ligands directly along the dissociation pathways and could exert long-range allosteric effects. © 2018 Wiley Periodicals, Inc.

Inexpensive Method for Selecting Receptor Structures for Virtual Screening.

This article introduces a screening performance index (SPI) to help select from a number of experimental structures one or a few that are more likely to identify more actives among its top hits from virtual screening of a compound library. It achieved this by docking only known actives to the experimental structures without considering a large number of decoys to reduce computational costs. The SPI is calculated by using the docking energies of the actives to all the receptor structures. We evaluated the performance of the SPI by applying it to study eight protein systems: fatty acid binding protein adipocyte FABP4, serine/threonine-protein kinase BRAF, beta-1 adrenergic receptor ADRB1, TGF-beta receptor type I TGFR1, adenosylhomocysteinase SAHH, thyroid hormone receptor beta-1 THB, phospholipase A2 group IIA PA2GA, and cytochrome P450 3a4 CP3A4. We found that the SPI agreed with the results from other popular performance metrics such as Boltzmann-Enhanced Discrimination Receiver Operator Characteristics (BEDROC), Robust Initial Enhancement (RIE), Area Under Accumulation Curve (AUAC), and Enrichment Factor (EF) but is less expensive to calculate. SPI also performed better than the best docking energy, the molecular volume of the bound ligand, and the resolution of crystal structure in selecting good receptor structures for virtual screening. The implications of these findings were further discussed in the context of ensemble docking, in situations when no experimental structure for the targeted protein was available, or under circumstances when quick choices of receptor structures need to be made before quantitative indexes such as the SPI and BEDROC can be calculated.

SRmapper: a fast and sensitive genome-hashing alignment tool.

UNLABELLED: Modern sequencing instruments have the capability to produce millions of short reads every day. The large number of reads produced in conjunction with variations between reads and reference genomic sequences caused both by legitimate differences, such as single-nucleotide polymorphisms and insertions/deletions (indels), and by sequencer errors make alignment a difficult and computationally expensive task, and many reads cannot be aligned. Here, we introduce a new alignment tool, SRmapper, which in tests using real data can align 10s of billions of base pairs from short reads to the human genome per computer processor day. SRmapper tolerates a higher number of mismatches than current programs based on Burrows-Wheeler transform and finds about the same number of alignments in 2-8× less time depending on read length (with higher performance gain for longer read length). The current version of SRmapper aligns both single and pair-end reads in base space fastq format and outputs alignments in Sequence Alignment/Map format. SRmapper uses a probabilistic approach to set a default number of mismatches allowed and determines alignment quality. SRmapper's memory footprint (∼2.5 GB) is small enough that it can be run on a computer with 4 GB of random access memory for a genome the size of a human. Finally, SRmapper is designed so that its function can be extended to finding small indels as well as long deletions and chromosomal translocations in future versions. AVAILABILITY: http://www.umsl.edu/∼wongch/software.html.

A new class of salicylic acid derivatives for inhibiting YopH of Yersinia pestis.

Previously, we identified a class of salicylic acid derivatives that display inhibitory activity against the protein tyrosine phosphatase YopH from Yersinia pestis. Because docking study suggested that the large phenyl ring attaching to the salicylic acid core might be exposed to the solvent and might not contribute significantly to binding, we have developed a new class of compounds that no longer contain this phenyl ring. We first devised a synthetic scheme for the compounds and then developed an automated computational screening model surrounding this synthetic scheme to help select a small number of compounds for synthesis and experimental testing. Based on this computational screening model and the analysis of the structure-activity relationship of our previous class of compounds, we have synthesized eight compounds and found five that yield micromolar activity. When applying in a larger scale, the synthetic scheme and the computational screening model developed here should help to identify even more potent inhibitors in the future.

Quantitative Prediction of Dissociation Rates of PYK2 Ligands Using Umbrella Sampling and Milestoning.

We used umbrella sampling and the milestoning simulation method to study the dissociation of multiple ligands from protein kinase PYK2. The activation barriers obtained from the potential of mean force of the umbrella sampling simulations correlated well with the experimental dissociation rates. Using the zero-temperature string method, we obtained optimized paths along the free-energy surfaces for milestoning simulations of three ligands with a similar structure. The milestoning simulations gave an absolute dissociation rate within 2 orders of magnitude of the experimental value for two ligands but at least 3 orders of magnitude too high for the third. Despite the similarity in their structures, the ligands took different pathways to exit from the binding site of PYK2, making contact with different sets of residues. In addition, the protein experienced different conformational changes for dissociation of the three ligands.

15 Years of molecular simulation of drug-binding kinetics.

INTRODUCTION: Drug-binding kinetics has been increasingly recognized as an important factor to be considered in drug discovery. Long residence time could prolong the action of some drugs while produce toxicity on others. Early evaluation of the binding kinetics of drug candidates could reduce attrition rate late in the drug discovery process. Computational prediction of drug-binding kinetics is useful as compounds can be evaluated even before they are made. However, simulation of drug-binding kinetics is a challenging problem because of the long-time scale involved. Nevertheless, significant progress has been made. AREAS COVERED: This review illustrates the rapid evolution of qualitative to quantitative molecular dynamics-based methods that have been developed over the last 15 years. EXPERT OPINION: The development of new methods based on molecular dynamics simulations now enables computation of absolute association/dissociation rate constants. Cheaper methods capable of identifying candidates with fast or slow binding kinetics, or rank-ordering rate constants are also available. Together, these methods have generated useful insights into the molecular mechanisms of drug-binding kinetics, and the design of drug candidates with therapeutically favorable kinetics. Although predicting absolute rate constants is still expensive and challenging, rapid improvement is expected in the coming years with the continuing refinement of current technologies, development of new methodologies, and the utilization of machine learning.

Simulation reveals two major docking pathways between the hexapeptide GDYMNM and the catalytic domain of the insulin receptor protein kinase.

Extending a previous mining-minima approach to identifying the docking pathways between the hexapeptide GDYMNM and the catalytic domain of the insulin receptor tyrosine kinase (IRK), we found two major docking pathways connecting the binding pocket and the surface of the protein. One pathway was more likely to lead to phosphate transfer from ATP to the peptide as the distance between the γ-phosphate of ATP and the hydroxyl oxygen of the target tyrosine approached one that could facilitate reaction. The movement of the peptide along the pathways was found to couple with residues in the activation loop of the protein. Although these residues might not affect binding affinity, they might influence the kinetics of peptide entrance and release.

Influence of kinetics of drug binding on EGFR signaling: a comparative study of three EGFR signaling pathway models.

We used three models of the epidermal growth factor receptor (EGFR) signaling pathway mimicking three different cell lines to study the effects of kinetics of drug binding on influencing molecular signaling in the pathways. With no incubation of drugs before the external cue epidermal growth factor (EGF) was applied, we found that fast kinetics of binding to protein kinases was advantageous in suppressing the production of the Extracellular signal-regulated kinase (ERK) that triggers cell proliferation, with some exceptions. Incubation of a drug with a protein kinase target for an hour before a pathway was initiated with an external cue made kinetics less significant, so did high concentration of drugs. In addition, we found that applying a drug to a protein kinase mostly affected downstream signaling although upstream events were also affected in a few cases. In examining whether applying two drugs to two protein kinase targets in the pathways could produce synergistic effects, we found positive, negative, or no effects, depending on the protein kinases targeted and the pathway model considered.

Qualitative Prediction of Ligand Dissociation Kinetics from Focal Adhesion Kinase Using Steered Molecular Dynamics.

Most early-stage drug discovery projects focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. Since many approved drugs have nonequilibrium binding characteristics, there has been increasing interest in optimizing binding kinetics early in the drug discovery process. As focal adhesion kinase (FAK) is an important drug target, we examine whether steered molecular dynamics (SMD) can be useful for identifying drug candidates with the desired drug-binding kinetics. In simulating the dissociation of 14 ligands from FAK, we find an empirical power-law relationship between the simulated time needed for ligand unbinding and the experimental rate constant for dissociation, with a strong correlation depending on the SMD force used. To improve predictions, we further develop regression models connecting experimental dissociation rate with various structural and energetic quantities derived from the simulations. These models can be used to predict dissociation rates from FAK for related compounds.

EDock-ML: A web server for using ensemble docking with machine learning to aid drug discovery.

EDock-ML is a web server that facilitates the use of ensemble docking with machine learning to help decide whether a compound is worthwhile to be considered further in a drug discovery process. Ensemble docking provides an economical way to account for receptor flexibility in molecular docking. Machine learning improves the use of the resulting docking scores to evaluate whether a compound is likely to be useful. EDock-ML takes a bottom-up approach in which machine-learning models are developed one protein at a time to improve predictions for the proteins included in its database. Because the machine-learning models are intended to be used without changing the docking and model parameters with which the models were trained, novice users can use it directly without worrying about what parameters to choose. A user simply submits a compound specified by an ID from the ZINC database (Sterling, T.; Irwin, J. J., J Chem Inf Model 2015, 55[11], 2,324-2,337.) or upload a file prepared by a chemical drawing program and receives an output helping the user decide the likelihood of the compound to be active or inactive for a drug target. EDock-ML can be accessed freely at edock-ml.umsl.edu.

Flexible ligand-flexible protein docking in protein kinase systems.

Quite a few reviews on molecular docking have already appeared. This mini-review focuses on methods that incorporate protein flexibility in docking rather than those that treat protein targets as rigid molecules. This is still a challenging problem but there are encouraging recent advances. These methods will be reviewed particularly in light of their applications to protein kinases and phosphatases. In addition to obtaining correct docking pose, recent developments on exploring docking pathways are also highlighted.

Conformational selection of protein kinase A revealed by flexible-ligand flexible-protein docking.

Protein kinases have high structural plasticity: their structure can change significantly, depending on what ligands are bound to them. Rigid-protein docking methods are not capable of describing such effects. Here, we present a new flexible-ligand flexible-protein docking model in which the protein can adopt conformations between two extremes observed experimentally. The model utilized a molecular dynamics-based simulated annealing cycling protocol and a distance-dependent dielectric model to perform docking. By testing this model on docking four diverse ligands to protein kinase A, we found that the ligands were able to dock successfully to the protein with the proper conformations of the protein induced. By imposing relatively soft conformational restraints to the protein during docking, this model reduced computational costs yet permitted essential conformational changes that were essential for these inhibitors to dock properly to the protein. For example, without adequate movement of the glycine-rich loop, it was difficult for the ligands to move from the surface of the protein to the binding site. In addition, these simulations called for better ways to compare simulation results with experiment other than using the popular root-mean-square deviation between the structure of a ligand in a docking pose and that in experiment because the structure of the protein also changed. In this work, we also calculated the correlation coefficient between protein-ligand/protein-protein distances in the docking structure and those in the crystal structure to check how well a ligand docked into the binding site of the protein and whether the proper conformation of the protein was induced.

A computational study of the phosphorylation mechanism of the insulin receptor tyrosine kinase.

Although various groups have studied the phosphorylation mechanism of the insulin receptor tyrosine kinase (IRK), an unanimous picture has not yet emerged. In this work, we performed a computational study to gain further insights. We first built a structural model of the reactant complex with the guide of several crystal structures and previous computational studies of the cyclic AMP-dependent protein kinase. We then optimized the structure by performing geometry optimization using a quantum mechanical model containing nearly 300 atoms. A reaction path was then traced between the reactant and the product by using a multiple coordinate-driven method. The calculations mapped out a sequence of structural changes depicting the conversion of the reactant to the product. Analysis of the structural changes revealed the formation of a dissociative transition state and the involvement of a proton transfer from the hydroxyl group of the tyrosyl residue of the peptide substrate to a conserved aspartate in the active site of the enzyme. The proton transfer began well before the transition state was reached and finished only shortly before the product was completely formed. In addition, the formation of a hydrogen bonding network among Arg1136, Asp1132, the gamma-phosphate of ATP, and the tyrosine residue of the substrate appeared to hold the latter two in a near-attack position for reaction. The model estimated a reaction barrier of 14 kcal/mol, semiquantitatively in accord with experiment.

Flexible protein-flexible ligand docking with disrupted velocity simulated annealing.

By docking flexible balanol to a rigid model of protein kinase A (PKA), we found that a new simulated annealing protocol termed disrupted velocity simulated annealing (DIVE-SA) outperformed the replica-exchange method and the traditional simulated annealing method in identifying the correct docking pose. In this protocol, the atomic velocities were reassigned periodically to encourage the system to sample a large conformational space. We also found that scaling potential energy surface to reduce structural transition barriers could further facilitate docking. The DIVE-SA method was then evaluated on its ability to perform flexible ligand-flexible protein docking of three ligands (balanol, a balanol analog, and ATP) to PKA. To reduce computational time and to avoid possible unphysical structural changes resulting from the use of nonoptimal force fields, a soft restrain was applied to keep the root-mean-square-deviation (RMSD) between instantaneous protein structures and a chosen reference structure small. Because the restrain was applied to the overall RMSD rather than to individual atoms, a protein could still experience relatively large conformational changes during docking. To examine the impact of applying such a restrain on docking, we constructed two semi-flexible protein models by choosing two different crystal structures as reference. Both the balanol analog and ATP were able to dock to either one of these semi-flexible protein models. On the other hand, balanol could only dock well to one of them. Further analysis indicated that the restrain on the glycine-rich loop was too strong, preventing it to adjust its structure to accommodate balanol in the binding pocket of PKA. Removing the restrain on the glycine-rich loop resulted in much better docking poses. This finding demonstrates the important role that the flexibility of the glycine-rich loop play in accepting different ligands and should profitably not be restrained in molecular docking so that more diverse ligands can be studied.

Variable atomic radii for continuum-solvent electrostatics calculation.

We have developed a method to improve the description of solute cavity defined by the interlocking-sphere model for continuum-solvent electrostatics calculations. Many models choose atomic radii from a finite set of atom types or uses an even smaller set developed by Bondi [J. Phys. Chem. 68, 441 (1964)]. The new model presented here allowed each atom to adapt its radius according to its chemical environment. This was achieved by first approximating the electron density of a molecule by a superposition of atom-centered spherical Gaussian functions. The parameters of the Gaussian functions were then determined by optimizing a function that minimized the difference between the properties from the model and those from ab initio quantum calculations. These properties included the electrostatics potential on molecular surface and the electron density within the core of each atom. The size of each atom was then determined by finding the radius at which the electron density associated with the atom fell to a prechosen value. This value was different for different chemical elements and was chosen such that the averaged radius for each chemical element in a training set of molecules matched its Bondi radius. Thus, our model utilized only a few adjustable parameters-the above density cutoff values for different chemical elements-but had the flexibility of allowing every atom to adapt its radius according to its chemical environment. This variable-radii model gave better solvation energy for 31 small neutral molecules than the Bondi radii did, especially for a quantum mechanics/Poisson-Boltzmann approach we developed earlier. The improvement was most significant for molecules with large dipole moment. Future directions for further improvement are also discussed.

A mining minima approach to exploring the docking pathways of p-nitrocatechol sulfate to YopH.

Using the docking of p-nitrocatechol sulfate to Yersinia protein tyrosine phosphatase YopH as an example, we showed that an approach based on mining minima followed by cluster and similarity analysis could generate useful insights into docking pathways. Our simulation treated both the ligand and the protein as flexible molecules so that the coupling between their motion could be properly accounted for. Our simulation identified three docking poses; the one with the lowest energy agreed well with experimental structure. The model also predicted the side-chain conformations of the amino acids lying in the binding pocket correctly with the exception of three residues that appeared to be stabilized by two structural water molecules in the crystal structure. The implicit solvent model employed in the simulation could not capture such effects well. We also found four major pathways leading to these docking poses after the ligand entered the mouth of the binding pocket. In addition, the sulfate group of p-nitrocatechol sulfate was found to be important both in binding the ligand to the pocket and in guiding the ligand to dock into the pocket. The coupling of the motion between the protein and the ligand also played an important role in facilitating ligand loading and unloading.