MixMD Probeview: Robust Binding Site Prediction from Cosolvent Simulations

Exploring Aromatic Cage Flexibility Using Cosolvent Molecular Dynamics Simulations─An In-Silico Case Study of Tudor Domains

Cosolvent molecular dynamics (MD) simulations have proven to be powerful in silico tools to predict hotspots for binding regions on protein surfaces. In the current study, the method was adapted and applied to two Tudor domain-containing proteins, namely Spindlin1 (SPIN1) and survival motor neuron protein (SMN). Tudor domains are characterized by so-called aromatic cages that recognize methylated lysine residues of protein targets. In the study, the conformational transitions from closed to open aromatic cage conformations were investigated by performing MD simulations with cosolvents using six different probe molecules. It is shown that a trajectory clustering approach in combination with volume and atomic distance tracking allows a reasonable discrimination between open and closed aromatic cage conformations and the docking of inhibitors yields very good reproducibility with crystal structures. Cosolvent MDs are suitable to capture the flexibility of aromatic cages and thus represent a promising tool for the optimization of inhibitors.

Unraveling the Ligand-Binding Sites of CYP3A4 by Molecular Dynamics Simulations with Solvent Probes

Cytochrome P450 3A4 (CYP3A4) is one of the most important drug-metabolizing enzymes in the human body and is well known for its complicated, atypical kinetic characteristics. The existence of multiple ligand-binding sites in CYP3A4 has been widely recognized as being capable of interfering with the active pocket through allosteric effects. The identification of ligand-binding sites other than the canonical active site above the heme is especially important for understanding the atypical kinetic characteristics of CYP3A4 and the intriguing association between the ligand and the receptor. In this study, we first employed mixed-solvent molecular dynamics (MixMD) simulations coupled with the online computational predictive tools to explore potential ligand-binding sites in CYP3A4. The MixMD approach demonstrates better performance in dealing with the receptor flexibility compared with other computational tools. From the sites identified by MixMD, we then picked out multiple sites for further exploration using ensemble docking and conventional molecular dynamics (cMD) simulations. Our results indicate that three extra sites are suitable for ligand binding in CYP3A4, including one experimentally confirmed site and two novel sites.

Colabind: A Cloud-Based Approach for Prediction of Binding Sites Using Coarse-Grained Simulations with Molecular Probes

Binding site prediction is a crucial step in understanding protein-ligand and protein-protein interactions (PPIs) with broad implications in drug discovery and bioinformatics. This study introduces Colabind, a robust, versatile, and user-friendly cloud-based approach that employs coarse-grained molecular dynamics simulations in the presence of molecular probes, mimicking fragments of drug-like compounds. Our method has demonstrated high effectiveness when validated across a diverse range of biological targets spanning various protein classes, successfully identifying orthosteric binding sites, as well as known druggable allosteric or PPI sites, in both experimentally determined and AI-predicted protein structures, consistently placing them among the top-ranked sites. Furthermore, we suggest that careful inspection of the identified regions with a high affinity for specific probes can provide valuable insights for the development of pharmacophore hypotheses. The approach is available at https://github.com/porekhov/CG_probeMD.

Applications of Molecular Dynamics Simulations in Drug Discovery

In the current drug development process, molecular dynamics (MD) simulations have proven to be very useful. This chapter provides an overview of the current applications of MD simulations in drug discovery, from detecting protein druggable sites and validating drug docking outcomes to exploring protein conformations and investigating the influence of mutations on its structure and functions. In addition, this chapter emphasizes various strategies to improve the conformational sampling efficiency in molecular dynamics simulations. With a growing computer power and developments in the production of force fields and MD techniques, the importance of MD simulations in helping the drug development process is projected to rise significantly in the future.

Application of Mixed-Solvent Molecular Dynamics Simulations for Prediction of Allosteric Sites on G Protein-Coupled Receptors

The development of small molecule allosteric modulators acting at G protein-coupled receptors (GPCRs) is becoming increasingly attractive. Such compounds have advantages over traditional drugs acting at orthosteric sites on these receptors, in particular target specificity. However, the number and locations of druggable allosteric sites within most clinically relevant GPCRs are unknown. In the present study, we describe the development and application of a mixed-solvent molecular dynamics (MixMD)-based method for the identification of allosteric sites on GPCRs. The method employs small organic probes with druglike qualities to identify druggable hotspots in multiple replicate short-timescale simulations. As proof of principle, we first applied the method retrospectively to a test set of five GPCRs (cannabinoid receptor type 1, C-C chemokine receptor type 2, M2 muscarinic receptor, P2Y purinoceptor 1, and protease-activated receptor 2) with known allosteric sites in diverse locations. This resulted in the identification of the known allosteric sites on these receptors. We then applied the method to the μ-opioid receptor. Several allosteric modulators for this receptor are known, although the binding sites for these modulators are not known. The MixMD-based method revealed several potential allosteric sites on the mu-opioid receptor. Implementation of the MixMD-based method should aid future efforts in the structure-based drug design of drugs targeting allosteric sites on GPCRs. SIGNIFICANCE STATEMENT: Allosteric modulation of G protein-coupled receptors (GPCRs) has the potential to provide more selective drugs. However, there are limited structures of GPCRs bound to allosteric modulators, and obtaining such structures is problematic. Current computational methods utilize static structures and therefore may not identify hidden or cryptic sites. Here we describe the use of small organic probes and molecular dynamics to identify druggable allosteric hotspots on GPCRs. The results reinforce the importance of protein dynamics in allosteric site identification.

Is the Stalk of the SARS-CoV-2 Spike Protein Druggable?

The spike protein is key to SARS-CoV-2 high infectivity because it facilitates the receptor binding domain (RBD) encounter with ACE2. As targeting subunit S1 has not yet delivered an ACE2-binding inhibitor, we have assessed the druggability of the conserved segment of the spike protein stalk within subunit S2 by means of an integrated computational approach that combines the molecular docking of an optimized library of fragments with high-throughput molecular dynamics simulations. The high propensity of the spike protein to mutate in key regions that are responsible for the recognition of the human angiotensin-converting enzyme 2 (hACE2) or for the recognition of antibodies, has made subunit S1 of the spike protein difficult to target. Despite the inherent flexibility of the stalk region, our results suggest two hidden interhelical binding sites, whose accessibility is only partially hampered by glycan residues.

Discovery of the Cryptic Sites of SARS-CoV-2 Papain-like Protease and Analysis of Its Druggability

In late 2019, a new coronavirus (CoV) caused the outbreak of a deadly respiratory disease, resulting in the COVID-19 pandemic. In view of the ongoing pandemic, there is an immediate need to find drugs to treat patients. SARS-CoV-2 papain-like cysteine protease (PLpro) not only plays an important role in the pathogenesis of the virus but is also a target protein for the development of inhibitor drugs. Therefore, to develop targeted inhibitors, it is necessary to analyse and verify PLpro sites and explore whether there are other cryptic binding pockets with better activity. In this study, first, we detected the site of the whole PLpro protein by sitemap of Schrödinger (version 2018), the cavity of LigBuilder V3, and DeepSite, and roughly judged the possible activated binding site area. Then, we used the mixed solvent dynamics simulation (MixMD) of probe molecules to induce conformational changes in the protein to find the possible cryptic active sites. Finally, the TRAPP method was used to predict the druggability of cryptic pockets and analyse the changes in the physicochemical properties of residues around these sites. This work will help promote the research of SARS-CoV-2 PLpro inhibitors.

Exploring the druggability of oxidized low-density lipoprotein (ox-LDL) receptor, LOX-1, a proatherogenic drug target involved in atherosclerosis

Lectin-like oxidized low-density lipoprotein (ox-LDL) receptor 1 (LOX-1) is a vital scavenger receptor involved in ox-LDL binding, internalization, and subsequent proatherogenic signaling leading to cellular dysfunction and atherosclerotic plaque formation. Existing data suggest that modulation of ox-LDL - LOX-1 interaction can prevent or slow down atherosclerosis. Therefore, we utilized computational methods such as multi-solvent simulation and characterized two top-ranked druggable sites. Using systematic molecular docking followed by atomistic molecular dynamics simulation, we have identified and shortlisted small molecules from the NCI library that target two key binding sites. We demonstrate, using surface plasmon resonance (SPR), that four of the shortlisted molecules bind one-on-one to the purified C-terminal domain (CTLD) of LOX-1 receptor with high affinity (K_D), ranging from 4.9 nM to 20.1 μM. Further, we performed WaterMap analysis to understand the role of individual water molecules in small molecule binding and the LOX-1-ligand complex stability. Our data clearly show that LOX-1 is druggable with small molecules. Our study provides strategies to identify novel inhibitors to attenuate ox-LDL - LOX-1 interaction.

Computational Identification of Possible Allosteric Sites and Modulators of the SARS-CoV-2 Main Protease

In this study, we target the main protease (M^pro) of the SARS-CoV-2 virus as it is a crucial enzyme for viral replication. Herein, we report three plausible allosteric sites on M^pro that can expand structure-based drug discovery efforts for new M^pro inhibitors. To find these sites, we used mixed-solvent molecular dynamics (MixMD) simulations, an efficient computational protocol that finds binding hotspots through mapping the surface of unbound proteins with 5% cosolvents in water. We have used normal mode analysis to support our claim of allosteric control for these sites. Further, we have performed virtual screening against the sites with 361 hits from M^pro screenings available through the National Center for Advancing Translational Sciences (NCATS). We have identified the NCATS inhibitors that bind to the remote sites better than the active site of M^pro, and we propose these molecules may be allosteric regulators of the system. After identifying our sites, new X-ray crystal structures were released that show fragment molecules in the sites we found, supporting the notion that these sites are accurate and druggable.

Application of Site-Identification by Ligand Competitive Saturation in Computer-Aided Drug Design

Site Identification by Ligand Competitive Saturation (SILCS) is a molecular simulation approach that uses diverse small solutes in aqueous solution to obtain functional group affinity patterns of a protein or other macromolecule. This involves employing a combined Grand Canonical Monte Carlo (GCMC)-molecular dynamics (MD) method to sample the full 3D space of the protein, including deep binding pockets and interior cavities from which functional group free energy maps (FragMaps) are obtained. The information content in the maps, which include contributions from protein flexibilty and both protein and functional group desolvation contributions, can be used in many aspects of the drug discovery process. These include identification of novel ligand binding pockets, including allosteric sites, pharmacophore modeling, prediction of relative protein-ligand binding affinities for database screening and lead optimization efforts, evaluation of protein-protein interactions as well as in the formulation of biologics-based drugs including monoclonal antibodies. The present article summarizes the various tools developed in the context of the SILCS methodology and their utility in computer-aided drug design (CADD) applications, showing how the SILCS toolset can improve the drug-development process on a number of fronts with respect to both accuracy and throughput representing a new avenue of CADD applications.

Divide and Conquer. Pocket-Opening Mixed-Solvent Simulations in the Perspective of Docking Virtual Screening Applications for Drug Discovery

The existence of a druggable binding pocket is a prerequisite for computational drug-target interaction studies including virtual screening. Retrospective studies have shown that extended sampling methods like Markov State Modeling and mixed-solvent simulations can identify cryptic pockets relevant for drug discovery. Here, we apply a combination of mixed-solvent molecular dynamics (MD) and time-structure independent component analysis (TICA) to four retrospective case studies: NPC2, the CECR2 bromodomain, TEM-1, and MCL-1. We compare previous experimental and computational findings to our results. It is shown that the successful identification of cryptic pockets depends on the system and the cosolvent probes. We used alternative TICA internal features such as the unbiased backbone coordinates or backbone dihedrals versus biased interatomic distances. We found that in the case of NPC2, TEM-1, and MCL-1, the use of unbiased features is able to identify cryptic pockets, although in the case of the CECR2 bromodomain, more specific features are required to properly capture a pocket opening. In the perspective of virtual screening applications, it is shown how docking studies with the parent ligands depend critically on the conformational state of the targets.

Fragment Hotspot Mapping to Identify Selectivity-Determining Regions between Related Proteins

Selectivity is a crucial property in small molecule development. Binding site comparisons within a protein family are a key piece of information when aiming to modulate the selectivity profile of a compound. Binding site differences can be exploited to confer selectivity for a specific target, while shared areas can provide insights into polypharmacology. As the quantity of structural data grows, automated methods are needed to process, summarize, and present these data to users. We present a computational method that provides quantitative and data-driven summaries of the available binding site information from an ensemble of structures of the same protein. The resulting ensemble maps identify the key interactions important for ligand binding in the ensemble. The comparison of ensemble maps of related proteins enables the identification of selectivity-determining regions within a protein family. We applied the method to three examples from the well-researched human bromodomain and kinase families, demonstrating that the method is able to identify selectivity-determining regions that have been used to introduce selectivity in past drug discovery campaigns. We then illustrate how the resulting maps can be used to automate comparisons across a target protein family.

Polyester Dendrimers Based on Bis-MPA for Doxorubicin Delivery

Although doxorubicin (DOX) is one of the most used chemotherapeutic drugs due to its efficacy against a wide group of cancer types, it presents severe side effects. As such, intensive research is being carried out to find new nanoscale systems that can help to overcome this problem. Polyester dendrimers based on the monomer 2,2-bis(hydroxymethyl)propionic acid (bis-MPA) are very promising systems for biomedical applications due to their biodegradability properties. In this study, bis-MPA-based dendrimers were, for the first time, evaluated as DOX delivery vehicles. Generations 4 and 5 of bis-MPA-based dendrimers with hydroxyl groups at the surface were used (B-G4-OH and B-G5-OH), together with dendrimers partially functionalized with amine groups (B-G4-NH₂/OH and B-G5-NH₂/OH). Partial functionalization was chosen because the main purpose was to compare the effect of different functional groups on dendrimers' drug delivery behavior without compromising cell viability, which is often affected by dendrimers' cationic charge. Results revealed that bis-MPA-based dendrimers were cytocompatible, independently of the chemical groups that were present at their surface. The B-G4-NH₂/OH and B-G5-NH₂/OH dendrimers were able to retain a higher number of DOX molecules, but the in vitro release of the drug was faster. On the contrary, the hydroxyl-terminated dendrimers exhibited a lower loading capacity but were able to deliver the drug in a more sustained manner. These results were in accordance with the cytotoxicity studies performed in several models of cancer cell lines and human mesenchymal stem cells. Overall, the results confirmed that it is possible to tune the drug delivery properties of bis-MPA-based dendrimers by modifying surface functionalization. Moreover, molecular modeling studies provided insights into the nature of the interactions established between the drug and the bis-MPA-based dendrimers─DOX molecules attach to their surface rather than being physically encapsulated.

Mixed-solvent molecular dynamics simulation-based discovery of a putative allosteric site on regulator of G protein signaling 4

Regulator of G protein signaling 4 (RGS4) is an intracellular protein that binds to the G_α subunit ofheterotrimeric G proteins and aids in terminating G protein coupled receptor signaling. RGS4 has been implicated in pain, schizophrenia, and the control of cardiac contractility. Inhibitors of RGS4 have been developed but bind covalently to cysteine residues on the protein. Therefore, we sought to identify alternative druggable sites on RGS4 using mixed-solvent molecular dynamics simulations, which employ low concentrations of organic probes to identify druggable hotspots on the protein. Pseudo-ligands were placed in consensus hotspots, and perturbation with normal mode analysis led to the identification and characterization of a putative allosteric site, which would be invaluable for structure-based drug design of non-covalent, small molecule inhibitors. Future studies on the mechanism of this allostery will aid in the development of novel therapeutics targeting RGS4.

Detection of Binding Sites on SARS-CoV-2 Spike Protein Receptor-Binding Domain by Molecular Dynamics Simulations in Mixed Solvents

The novel SARS-CoV-2 uses ACE2 (Angiotensin-Converting Enzyme 2) receptor as an entry point. Insights on S protein receptor-binding domain (RBD) interaction with ACE2 receptor and drug repurposing has accelerated drug discovery for the novel SARS-CoV-2 infection. Finding small molecule binding sites in S protein and ACE2 interface is crucial in search of effective drugs to prevent viral entry. In this study, we employed molecular dynamics simulations in mixed solvents together with virtual screening to identify small molecules that could be potential inhibitors of S protein -ACE2 interaction. Observation of organic probe molecule localization during the simulations revealed multiple sites at the S protein surface related to small molecule, antibody, and ACE2 binding. In addition, a novel conformation of the S protein was discovered that could be stabilized by small molecules to inhibit attachment to ACE2. The most promising binding site on RBD-ACE2 interface was targeted with virtual screening and top-ranked compounds (DB08248, DB02651, DB03714, and DB14826) are suggested for experimental testing. The protocol described here offers an extremely fast method for characterizing key proteins of a novel pathogen and for the identification of compounds that could inhibit or accelerate spreading of the disease.

Finding Druggable Sites in Proteins Using TACTICS

Structure-based drug discovery efforts require knowledge of where drug-binding sites are located on target proteins. To address the challenge of finding druggable sites, we developed a machine-learning algorithm called TACTICS (trajectory-based analysis of conformations to identify cryptic sites), which uses an ensemble of molecular structures (such as molecular dynamics simulation data) as input. First, TACTICS uses k-means clustering to select a small number of conformations that represent the overall conformational heterogeneity of the data. Then, TACTICS uses a random forest model to identify potentially bindable residues in each selected conformation, based on protein motion and geometry. Lastly, residues in possible binding pockets are scored using fragment docking. As proof-of-principle, TACTICS was applied to the analysis of simulations of the SARS-CoV-2 main protease and methyltransferase and the Yersinia pestis aryl carrier protein. Our approach recapitulates known small-molecule binding sites and predicts the locations of sites not previously observed in experimentally determined structures. The TACTICS code is available at https://github.com/Albert-Lau-Lab/tactics_protein_analysis.

EXPRORER: Rational Cosolvent Set Construction Method for Cosolvent Molecular Dynamics Using Large-Scale Computation

Cosolvent molecular dynamics (CMD) simulations involve an MD simulation of a protein in the presence of explicit water molecules mixed with cosolvent molecules to perform hotspot detection, binding site identification, and binding energy estimation, while other existing methods (e.g., MixMD, SILCS, and MDmix) utilize small molecules that represent functional groups of compounds. However, the cosolvent selections employed in these methods differ and there are only a few cosolvents that are commonly used in these methods. In this study, we proposed a systematic method for constructing a set of cosolvents for drug discovery, termed the EXtended PRObes set construction by REpresentative Retrieval (EXPRORER). First, we extracted typical substructures from FDA-approved drugs, generated 138 cosolvent structures, and for each cosolvent molecule, we conducted CMD simulations to generate a spatial probability distribution map of cosolvent atoms (PMAP). Analyses of PMAP similarity revealed that a cosolvent pair with a PMAP similarity greater than 0.70-0.75 shared similar structural features. We present a method for the construction of a cosolvent subset that satisfies a similarity threshold for all cosolvents, and we tested the constructed sets for four proteins. To our knowledge, this is the first study to include a systematic proposal for cosolvent set construction, and thus, the EXPRORER cosolvents will provide deeper insights into ligand binding sites of various proteins.

Analysis of the Binding Sites on BAX and the Mechanism of BAX Activators through Extensive Molecular Dynamics Simulations

The BAX protein is a pro-apoptotic member of the Bcl-2 family, which triggers apoptosis by causing permeabilization of the mitochondrial outer membrane. However, the activation mechanism of BAX is far from being understood. Although a few small-molecule BAX activators have been reported in the literature, their crystal structures in complex with BAX have not been resolved. So far, their binding modes were modeled at most by simple molecular docking efforts. Lack of an in-depth understanding of the activation mechanism of BAX hinders the development of more effective BAX activators. In this work, we employed cosolvent molecular dynamics simulation to detect the potential binding sites on the surface of BAX and performed a long-time molecular dynamics simulation (50 μs in total) to derive the possible binding modes of three BAX activators (i.e., BAM7, BTC-8, and BTSA1) reported in the literature. Our results indicate that the trigger, S184, and vMIA sites are the three major binding sites on the full-length BAX structure. Moreover, the canonical hydrophobic groove is clearly detected on the α9-truncated BAX structure, which is consistent with the outcomes of relevant experimental studies. Interestingly, it is observed that solvent probes bind to the trigger bottom pocket more stably than the PPI trigger site. Each activator was subjected to unbiased molecular dynamics simulations started at the three major binding sites in five parallel jobs. Our MD results indicate that all three activators tend to stay at the trigger site with favorable MM-GB/SA binding energies. BAM7 and BTSA1 can enter the trigger bottom pocket and thereby enhance the movement of the α1-α2 loop, which may be a key factor at the early stage of BAX activation. Our molecular modeling results may provide useful guidance for designing smart biological experiments to further explore BAX activation and directing structure-based efforts toward discovering more effective BAX activators.

Simultaneous Identification of Multiple Binding Sites in Proteins: A Statistical Mechanics Approach

We present an extension of the Poisson-Boltzmann model in which the solute of interest is immersed in an assembly of self-orienting Langevin water dipoles, anions, cations, and hydrophobic molecules, all of variable densities. Interactions between charges are controlled by electrostatics, while hydrophobic interactions are modeled with a Yukawa potential. We impose steric constraints by assuming that the system is represented on a cubic lattice. We also assume incompressibility; i.e., all sites of the lattice are occupied. This model, which we refer to as the Hydrophobic Dipolar Poisson-Boltzmann Langevin (HDPBL) model, leads to a system of two equations whose solutions give the water dipole, salt, and hydrophobic molecule densities, all of them in the presence of the others in a self-consistent way. We use those to study the organization of the ions, cosolvent, and solvent molecules around proteins. In particular, peaks of densities are expected to reveal, simultaneously, the presence of compatible binding sites of different kinds on a protein. We have tested and validated the ability of HDPBL to detect pockets in proteins that bind to hydrophobic ligands, polar ligands, and charged small probes as well as to characterize the binding sites of lipids for membrane proteins.

Protein druggability assessment for natural products using in silico simulation: A case study with estrogen receptor and the flavonoid genistein

Traditional herbal medicine (THM) comprises a vast number of natural compounds. Most of them are metabolized into different structures after administration, which makes the clarification of THM's mode of action more complicated. To evaluate the biological activities of those components and metabolites, in silico simulation technology is helpful. We focused on mixed-solvent molecular dynamics (MD) simulation for druggability assessment of natural products. Mixed-solvent MD is an in silico simulation method for the exploration of ligand-binding sites on target proteins, which uses water and an organic molecule mixture. The selection of organic small molecules is an important factor for predicting the characteristics of natural products. In this study, we used the known crystal structure of estrogen receptors with genistein as a test case and explored fragments reflecting the characteristics of natural products. We found that structures with a 4-pyrone structure are more often included in the natural products database compared with the DrugBank database, and we selectively detected the known-binding sites of estrogen receptor α and β. The results indicate that the 4-pyrone structure might be promising for predicting the protein druggability of flavonoids. Additionally, mixed-solvent MD simulation discriminates the selectivity of genistein between estrogen receptor β and α, indicating that the simulation can be evaluated using indices that differ from those of traditional ligand docking. Although this approach is still in its early stages, it has the potential to provide valuable information for understanding the diverse biological activities of natural products.

Identification of Cryptic Binding Sites Using MixMD with Standard and Accelerated Molecular Dynamics

Protein dynamics play an important role in small molecule binding and can pose a significant challenge in the identification of potential binding sites. Cryptic binding sites have been defined as sites which require significant rearrangement of the protein structure to become physically accessible to a ligand. Mixed-solvent MD (MixMD) is a computational protocol which maps the surface of the protein using molecular dynamics (MD) of the unbound protein solvated in a 5% box of probe molecules with explicit water. This method has successfully identified known active and allosteric sites which did not require reorganization. In this study, we apply the MixMD protocol to identify known cryptic sites of 12 proteins characterized by a wide range of conformational changes. Of these 12 proteins, three require reorganization of side chains, five require loop movements, and four require movement of more significant structures such as whole helices. In five cases, we find that standard MixMD simulations are able to map the cryptic binding sites with at least one probe type. In two cases (guanylate kinase and TIE-2), accelerated MD, which increases sampling of torsional angles, was necessary to achieve mapping of portions of the cryptic binding site missed by standard MixMD. For more complex systems where movement of a helix or domain is necessary, MixMD was unable to map the binding site even with accelerated dynamics, possibly due to the limited timescale (100 ns for individual simulations). In general, similar conformational dynamics are observed in water-only simulations and those with probe molecules. This could imply that the probes are not driving opening events but rather take advantage of mapping sites that spontaneously open as part of their inherent conformational behavior. Finally, we show that docking to an ensemble of conformations from the standard MixMD simulations performs better than docking the apo crystal structure in nine cases and even better than half of the bound crystal structures. Poorer performance was seen in docking to ensembles of conformations from the accelerated MixMD simulations.

Straightforward Incorporation of Multiple Ligand Types into Molecular Dynamics Simulations for Efficient Binding Site Detection and Characterization

Binding site identification and characterization is an important initial step in structure-based drug design. To account for the effects of protein flexibility and solvation, several cosolvent molecular dynamics (MD) simulation methods that incorporate small organic molecules into the protein's solvent box to probe for binding sites have been developed. However, most of these methods are highly inefficient, as they allow for the use of only one probe type at a time, which means that multiple sets of simulations have to be performed to map different types of binding sites. The high probe concentrations used in some of these methods also necessitate the use of artificial repulsive forces to prevent the probes from aggregating. Here, we present multiple-ligand-mapping MD (mLMMD), a method that incorporates multiple types of probes for simultaneous and efficient mapping of different types of binding sites without the need for introduction of artificial forces that may cause unintended mapping artifacts. We validate the method on a diverse set of 10 proteins and show that the mLMMD probes are able to reliably identify hydrophobic, hydrogen-bonding, charged, and cryptic binding sites in all of the test cases. Our results also highlight the potential utility of mLMMD for virtual screening and rational drug design.

Advances in Molecular Dynamics Simulations and Enhanced Sampling Methods for the Study of Protein Systems

Molecular dynamics (MD) simulation is a rigorous theoretical tool that when used efficiently could provide reliable answers to questions pertaining to the structure-function relationship of proteins. Data collated from protein dynamics can be translated into useful statistics that can be exploited to sieve thermodynamics and kinetics crucial for the elucidation of mechanisms responsible for the modulation of biological processes such as protein-ligand binding and protein-protein association. Continuous modernization of simulation tools enables accurate prediction and characterization of the aforementioned mechanisms and these qualities are highly beneficial for the expedition of drug development when effectively applied to structure-based drug design (SBDD). In this review, current all-atom MD simulation methods, with focus on enhanced sampling techniques, utilized to examine protein structure, dynamics, and functions are discussed. This review will pivot around computer calculations of protein-ligand and protein-protein systems with applications to SBDD. In addition, we will also be highlighting limitations faced by current simulation tools as well as the improvements that have been made to ameliorate their efficiency.

Impact of electronic polarizability on protein-functional group interactions

Interactions of proteins with functional groups are key to their biological functions, making it essential that they be accurately modeled. To investigate the impact of the inclusion of explicit treatment of electronic polarizability in force fields on protein-functional group interactions, the additive CHARMM and Drude polarizable force field are compared in the context of the Site-Identification by Ligand Competitive Saturation (SILCS) simulation methodology from which functional group interaction patterns with five proteins for which experimental binding affinities of multiple ligands are available, were obtained. The explicit treatment of polarizability produces significant differences in the functional group interactions in the ligand binding sites including overall enhanced binding of functional groups to the proteins. This is associated with variations of the dipole moments of solutes representative of functional groups in the binding sites relative to aqueous solution with higher dipole moments systematically occurring in the latter, though exceptions occur with positively charged methylammonium. Such variation indicates the complex, heterogeneous nature of the electronic environments of ligand binding sites and emphasizes the inherent limitation of fixed charged, additive force fields for modeling ligand-protein interactions. These effects yield more defined orientation of the functional groups in the binding pockets and a small, but systematic improvement in the ability of the SILCS method to predict the binding orientation and relative affinities of ligands to their target proteins. Overall, these results indicate that the physical model associated with the explicit treatment of polarizability along with the presence of lone pairs in a force field leads to changes in the nature of the interactions of functional groups with proteins versus that occurring with additive force fields, suggesting the utility of polarizable force fields in obtaining a more realistic understanding of protein-ligand interactions.

Monte Carlo simulations using PELE to identify a protein-protein inhibitor binding site and pose

In silico binding site location and pose prediction for a molecule targeted at a large protein surface is a challenging task. We report a blind test with two peptidomimetic molecules that bind the flu virus hemagglutinin (HA) surface antigen, JNJ7918 and JNJ4796 (recently disclosed in van Dongen et al., Science, 2019, 363). Tests with a series of conventional approaches such as rigid (receptor) docking against available X-ray crystal structures or against an ensemble of structures generated by quick methodologies (NMA, homology modeling) gave mixed results, due to the shallowness and flexibility of the binding site and the sheer size of the target. However, tests with our Monte Carlo platform PELE in two protocols involving either exploration of the whole protein surface (global exploration), or the latter followed by refinement of best solutions (local exploration) yielded remarkably good results by locating the actual binding site and generating binding modes that recovered all native contacts found in the X-ray structures. Thus, the Monte Carlo scheme of PELE seems promising as a quick methodology to overcome the challenge of identifying entirely unknown binding sites and modes for protein–protein disruptors.

PELE prospectively unveils the binding site and mode of a protein–protein disruptor.

Cosolvent Analysis Toolkit (CAT): a robust hotspot identification platform for cosolvent simulations of proteins to expand the druggable proteome

Cosolvent Molecular Dynamics (MD) simulations are increasingly popular techniques developed for prediction and characterization of allosteric and cryptic binding sites, which can be rendered “druggable” by small molecule ligands. Despite their conceptual simplicity and effectiveness, the analysis of cosolvent MD trajectories relies on pocket volume data, which requires a high level of manual investigation and may introduce a bias. In this work, we present CAT (Cosolvent Analysis Toolkit): an open-source, freely accessible analytical tool, suitable for automated analysis of cosolvent MD trajectories. CAT is compatible with commonly used molecular graphics software packages such as UCSF Chimera and VMD. Using a novel hybrid empirical force field scoring function, CAT accurately ranks the dynamic interactions between the macromolecular target and cosolvent molecules. To benchmark, CAT was used for three validated protein targets with allosteric and orthosteric binding sites, using five chemically distinct cosolvent molecules. For all systems, CAT has accurately identified all known sites. CAT can thus assist in computational studies aiming at identification of protein “hotspots” in a wide range of systems. As an easy-to-use computational tool, we expect that CAT will contribute to an increase in the size of the potentially ‘druggable’ human proteome.

Pharmmaker: Pharmacophore modeling and hit identification based on druggability simulations

Recent years have seen progress in druggability simulations, that is, molecular dynamics simulations of target proteins in solutions containing drug‐like probe molecules to characterize their drug‐binding abilities, if any. An important consecutive step is to analyze the trajectories to construct pharmacophore models (PMs) to use for virtual screening of libraries of small molecules. While considerable success has been observed in this type of computer‐aided drug discovery, a systematic tool encompassing multiple steps from druggability simulations to pharmacophore modeling, to identifying hits by virtual screening of libraries of compounds, has been lacking. We address this need here by developing a new tool, Pharmmaker, building on the DruGUI module of our ProDy application programming interface. Pharmmaker is composed of a suite of steps: (Step 1) identification of high affinity residues for each probe molecule type; (Step 2) selecting high affinity residues and hot spots in the vicinity of sites identified by DruGUI; (Step 3) ranking of the interactions between high affinity residues and specific probes; (Step 4) obtaining probe binding poses and corresponding protein conformations by collecting top‐ranked snapshots; and (Step 5) using those snapshots for constructing PMs. The PMs are then used as filters for identifying hits in structure‐based virtual screening. Pharmmaker, accessible online at http://prody.csb.pitt.edu/pharmmaker/, can be used in conjunction with other tools available in ProDy.

Cosolvent-Enhanced Sampling and Unbiased Identification of Cryptic Pockets Suitable for Structure-Based Drug Design

Modulating protein activity with small-molecules binding to cryptic pockets offers great opportunities to overcome hurdles in drug design. Cryptic sites are atypical binding sites in proteins that are closed in the absence of a stabilizing ligand and are thus inherently difficult to identify. Many studies have proposed methods to predict cryptic sites. However, a general approach to prospectively sample open conformations of these sites and to identify cryptic pockets in an unbiased manner suitable for structure-based drug design remains elusive. Here, we describe an all-atom, explicit cosolvent, molecular dynamics (MD) simulations-based workflow to sample the open states of cryptic sites and identify opened pockets, in a manner that does not require a priori knowledge about these sites. Furthermore, the workflow relies on a target-independent parametrization that only distinguishes between binding pockets for peptides or small molecules. We validated our approach on a diverse test set of seven proteins with crystallographically determined cryptic sites. The known cryptic sites were found among the three highest-ranked predicted cryptic sites, and an open site conformation was sampled and selected for most of the systems. Crystallographic ligand poses were well reproduced by docking into these identified open conformations for five of the systems. When the fully open state could not be reproduced, we were still able to predict the location of the cryptic site, or identify other cryptic sites that could be retrospectively validated with knowledge of the protein target. These characteristics render our approach valuable for investigating novel protein targets without any prior information.