Mixed-Probe Simulation and Probe-Derived Surface Topography Map Analysis for Ligand Binding Site Identification

Mechanisms of allostery at the viral surface through the eyes of molecular simulation

The outermost surface layer of any virus is formed by either a capsid shell or envelope. Such layers have traditionally been thought of as immovable structures, but it is becoming apparent that they cannot be viewed exclusively as static architectures protecting the viral genome. A limited number of proteins on the virion surface must perform a multitude of functions in order to orchestrate the viral life cycle, and allostery can regulate their structures at multiple levels of organization, spanning individual molecules, protomers, large oligomeric assemblies, or entire viral surfaces. Here, we review recent contributions from the molecular simulation field to viral surface allostery, with a particular focus on the trimeric spike glycoprotein emerging from the coronavirus surface, and the icosahedral flaviviral envelope complex. As emerging viral pathogens continue to pose a global threat, an improved understanding of viral dynamics and allosteric regulation will prove crucial in developing novel therapeutic strategies.

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.

Small-Molecule Inhibition of KRAS through Conformational Selection

Mutations in KRAS account for about 20% of human cancers. Despite the major progress in recent years toward the development of KRAS inhibitors, including the discovery of covalent inhibitors of the G12C KRAS variant for the treatment of non-small-cell lung cancer, much work remains to be done to discover broad-acting inhibitors to treat many other KRAS-driven cancers. In a previous report, we showed that a 308.4 Da small-molecule ligand [(2R)-2-(N'-(1H-indole-3-carbonyl)hydrazino)-2-phenyl-acetamide] binds to KRAS with low micro-molar affinity [Chem. Biol. Drug Des.2019; 94(2):1441-1456]. Binding of this ligand, which we call ACA22, to the p1 pocket of KRAS and its interactions with residues at beta-strand 1 and the switch loops have been supported by data from nuclear magnetic resonance spectroscopy and microscale thermophoresis experiments. However, the inhibitory potential of the compound was not demonstrated. Here, we show that ACA22 inhibits KRAS-mediated signal transduction in cells expressing wild type (WT) and G12D mutant KRAS and reduces levels of guanosine triphosphate-loaded WT KRAS more effectively than G12D KRAS. We ruled out the direct effect on nucleotide exchange or effector binding as possible mechanisms of inhibition using a variety of biophysical assays. Combining these observations with binding data that show comparable affinities of the compound for the active and inactive forms of the mutant but not the WT, we propose conformational selection as a possible mechanism of action of ACA22.

Uncovering cryptic pockets in the SARS-CoV-2 spike glycoprotein

The COVID-19 pandemic has prompted a rapid response in vaccine and drug development. Herein, we modeled a complete membrane-embedded SARS-CoV-2 spike glycoprotein and used molecular dynamics simulations with benzene probes designed to enhance discovery of cryptic pockets. This approach recapitulated lipid and host metabolite binding sites previously characterized by cryo-electron microscopy, revealing likely ligand entry routes, and uncovered a novel cryptic pocket with promising druggable properties located underneath the 617-628 loop. A full representation of glycan moieties was essential to accurately describe pocket dynamics. A multi-conformational behavior of the 617-628 loop in simulations was validated using hydrogen-deuterium exchange mass spectrometry experiments, supportive of opening and closing dynamics. The pocket is the site of multiple mutations associated with increased transmissibility found in SARS-CoV-2 variants of concern including Omicron. Collectively, this work highlights the utility of the benzene mapping approach in uncovering potential druggable sites on the surface of SARS-CoV-2 targets.

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.

Cosolvent and Dynamic Effects in Binding Pocket Search by Docking Simulations

The lack of conformational sampling in virtual screening projects can lead to inefficient results because many of the potential drugs may not be able to bind to the target protein during the static docking simulations. Here, we performed ensemble docking for around 2000 United States Food and Drug Administration (FDA)-approved drugs with the RNA-dependent RNA polymerase (RdRp) protein of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) as a target. The representative protein structures were generated by clustering classical molecular dynamics trajectories, which were evolved using three solvent scenarios, namely, pure water, benzene/water and phenol/water mixtures. The introduction of dynamic effects in the theoretical model showed improvement in docking results in terms of the number of strong binders and binding sites in the protein. Some of the discovered pockets were found only for the cosolvent simulations, where the nonpolar probes induced local conformational changes in the protein that lead to the opening of transient pockets. In addition, the selection of the ligands based on a combination of the binding free energy and binding free energy gap between the best two poses for each ligand provided more suitable binders than the selection of ligands based solely on one of the criteria. The application of cosolvent molecular dynamics to enhance the sampling of the configurational space is expected to improve the efficacy of virtual screening campaigns of future drug discovery projects.

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.

Approaches to inhibiting oncogenic K-Ras

Activating somatic K-Ras mutations are associated with >15% all human tumors and up to 90% of specific tumor types such as pancreatic cancer. Successfully inhibiting abnormal K-Ras signaling would therefore be a game changer in cancer therapy. However, K-Ras has long been considered an undruggable target for various reasons. This view is now changing by the discovery of allosteric inhibitors that directly target K-Ras and inhibit its functions, and by the identification of new mechanisms to dislodge it from the plasma membrane and thereby abrogate its cellular activities. In this review, we will discuss recent progresses and challenges to inhibiting aberrant K-Ras functions by these two approaches. We will also provide a broad overview of other approaches such as inhibition of K-Ras effectors, and offer a brief perspective on the way forward.

How to make an undruggable enzyme druggable: lessons from ras proteins

Significant advances have been made toward discovering allosteric inhibitors for challenging drug targets such as the Ras family of membrane-associated signaling proteins. Malfunction of Ras proteins due to somatic mutations is associated with up to a quarter of all human cancers. Computational techniques have played critical roles in identifying and characterizing allosteric ligand-binding sites on these proteins, and to screen ligand libraries against those sites. These efforts, combined with a wide range of biophysical, structural, biochemical and cell biological experiments, are beginning to yield promising inhibitors to treat malignancies associated with mutated Ras proteins. In this chapter, we discuss some of these developments and how the lessons learned from Ras might be applied to similar other challenging drug targets.

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.

A Benzene-Mapping Approach for Uncovering Cryptic Pockets in Membrane-Bound Proteins

Molecular dynamics (MD) simulations in combination with small organic probes present in the solvent have previously been used as a method to reveal cryptic pockets that may not have been identified in experimental structures. We report such a method implemented within the CHARMM force field using the GROMACS simulation package to effectively explore cryptic pockets on the surfaces of membrane-embedded proteins using benzene as a probe molecule. This method, for which we have made implementation files freely available, relies on modified nonbonded parameters in addition to repulsive potentials between membrane lipids and benzene molecules. The method was tested on part of the outer shell of the dengue virus (DENV), for which research into a safe and effective neutralizing antibody or drug molecule is still ongoing. In particular, the envelope (E) protein, associated with the membrane (M) protein, is a lipid membrane-embedded complex which forms a dimer in the mature viral envelope. Solvent mapping was performed for the full, membrane-embedded EM protein complex and compared with similar calculations performed for the isolated, soluble E protein ectodomain dimer in the solvent. Ectodomain-only simulations with benzene exhibited unfolding effects not observed in the more physiologically relevant membrane-associated systems. A cryptic pocket which has been experimentally shown to bind n-octyl-β-d-glucoside detergent was consistently revealed in all benzene-containing simulations. The addition of benzene also enhanced the flexibility and hydrophobic exposure of cryptic pockets at a key, functional interface in the E protein and revealed a novel, potentially druggable pocket that may be targeted to prevent conformational changes associated with viral entry into the cell.

Making NSCLC Crystal Clear: How Kinase Structures Revolutionized Lung Cancer Treatment

The parallel advances of different scientific fields provide a contemporary scenario where collaboration is not a differential, but actually a requirement. In this context, crystallography has had a major contribution on the medical sciences, providing a “face” for targets of diseases that previously were known solely by name or sequence. Worldwide, cancer still leads the number of annual deaths, with 9.6 million associated deaths, with a major contribution from lung cancer and its 1.7 million deaths. Since the relationship between cancer and kinases was unraveled, these proteins have been extensively explored and became associated with drugs that later attained blockbuster status. Crystallographic structures of kinases related to lung cancer and their developed and marketed drugs provided insight on their conformation in the absence or presence of small molecules. Notwithstanding, these structures were also of service once the initially highly successful drugs started to lose their effectiveness in the emergence of mutations. This review focuses on a subclassification of lung cancer, non-small cell lung cancer (NSCLC), and major oncogenic driver mutations in kinases, and how crystallographic structures can be used, not only to provide awareness of the function and inhibition of these mutations, but also how these structures can be used in further computational studies aiming at addressing these novel mutations in the field of personalized medicine.

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.

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.

Probing the Conformational and Energy Landscapes of KRAS Membrane Orientation

Membrane reorientation of oncogenic RAS proteins is emerging as an important modulator of their functions. Previous studies have shown that the most common orientations include those with either the three C-terminal α-helices (OS1) or N-terminal β-strands (OS2) of the catalytic domain facing the membrane. OS1 and OS2 differ by the degree to which the effector-interacting surface is occluded by the membrane. However, the relative stability of these states and the rates of transition between them remained undetermined. How mutations might modulate preferences for specific orientation states is also far from clear. The current work attempted to address these questions through a comprehensive analysis of two 20 μs-long atomistic molecular dynamics simulations. The simulations were conducted on the oncogenic G12D and Q61H KRAS mutants bound to an anionic lipid bilayer. G12D and Q61H are among the most prevalent cancer-causing mutations at the P-loop and switch 2 regions of KRAS, respectively. We found that both mutants fluctuate in a similar manner between OS1 and OS2 via an intermediate orientation OS0, and both favor the signaling competent OS1 and OS0 over the occluded OS2. However, they differ in the details, such as in the extent to which they sample OS1. Analysis of the orientation free-energy landscapes estimated from the simulations indicate that OS1 and OS2 are the most stable states. However, the overall free energy surface is rugged, indicating a large diversity of conformations including at least two substates in each orientation state that differ in stability only by about 0.5-1.0 kcal/mol. Reversible transitions between OS1 and OS2 occur via two well-defined pathways that traverse OS0. In the minimum energy path, helix 4 remains close to the membrane as the angle of the catalytic domain from the membrane plane changes, resulting in a barrier of ∼1 kcal/mol for OS1/OS2 interconversions. Estimation of the rates of the various transitions based on survival probabilities yielded two rate constants in the order of 10⁷ and 10⁶ s^-1, which we attribute to intrinsic protein conformational dynamics and transient protein-lipid interactions, respectively. The faster process dominates every transition, confirming a previous suggestion that RAS membrane reorientation is driven by conformational fluctuations rather than protein-lipid interactions.

Conformational and Dynamical Effects of Tyr32 Phosphorylation in K-Ras: Molecular Dynamics Simulation and Markov State Models Analysis

Phosphorylation of tyrosine 32 in K-Ras has been shown to influence function by disrupting the GTPase cycle. To shed light on the underlying mechanism and atomic basis of this process, we carried out a comparative investigation of the oncogenic G12D K-Ras mutant and its phosphorylated variant (pTyr32) using all-atom molecular dynamics simulations and Markov state models. We show that, despite sharing a number of common features, G12D and pTyr32-G12D K-Ras exhibit some distinct conformational states and fluctuations. In addition to notable differences in conformation and dynamics of residues surrounding the GTP binding site, nonlocal changes were observed at a number of loops. Switch I is more flexible in pTyr32-G12D K-Ras while switch II is more flexible in G12D K-Ras. We also used time-lagged independent component analysis and k-means clustering to identify five metastable states for each system. We utilized transition path theory to calculate the transition probabilities for each state to build a Markov state model for each system. These models and other close inspections suggest that the phosphorylation of Tyr32 strongly affects protein dynamics and the active site conformation, especially with regards to the canonical switch conformations and dynamics.

Solvents to Fragments to Drugs: MD Applications in Drug Design

Simulations of molecular dynamics (MD) are playing an increasingly important role in structure-based drug discovery (SBDD). Here we review the use of MD for proteins in aqueous solvation, organic/aqueous mixed solvents (MDmix) and with small ligands, to the classic SBDD problems: Binding mode and binding free energy predictions. The simulation of proteins in their condensed state reveals solvent structures and preferential interaction sites (hot spots) on the protein surface. The information provided by water and its cosolvents can be used very effectively to understand protein ligand recognition and to improve the predictive capability of well-established methods such as molecular docking. The application of MD simulations to the study of the association of proteins with drug-like compounds is currently only possible for specific cases, as it remains computationally very expensive and labor intensive. MDmix simulations on the other hand, can be used systematically to address some of the common tasks in SBDD. With the advent of new tools and faster computers we expect to see an increase in the application of mixed solvent MD simulations to a plethora of protein targets to identify new drug candidates.

Allosteric Modulation of Intact γ-Secretase Structural Dynamics

As a protease complex involved in the cleavage of amyloid precursor proteins that lead to the formation of amyloid β fibrils implicated in Alzheimer's disease, γ-secretase is an important target for developing therapeutics against Alzheimer's disease. γ-secretase is composed of four subunits: nicastrin (NCT) in the extracellular (EC) domain, presenilin-1 (PS1), anterior pharynx defective 1, and presenilin enhancer 2 in the transmembrane (TM) domain. NCT and PS1 play important roles in binding amyloid β precursor proteins and modulating PS1 catalytic activity. Yet, the molecular mechanisms of coupling between substrate/modulator binding and catalytic activity remain to be elucidated. Recent determination of intact human γ-secretase cryo-electron microscopy structure has opened the way for a detailed investigation of the structural dynamics of this complex. Our analysis, based on a membrane-coupled anisotropic network model, reveals two types of NCT motions, bending and twisting, with respect to PS1. These underlie the fluctuations between the "open" and "closed" states of the lid-like NCT with respect to a hydrophilic loop 1 (HL1) on PS1, thus allowing or blocking access of the substrate peptide (EC portion) to HL1 and to the neighboring helix TM2. In addition to this alternating access mechanism, fluctuations in the volume of the PS1 central cavity facilitate the exposure of the catalytic site for substrate cleavage. Druggability simulations show that γ-secretase presents several hot spots for either orthosteric or allosteric inhibition of catalytic activity, consistent with experimental data. In particular, a hinge region at the interface between the EC and TM domains, near the interlobe groove of NCT, emerges as an allo-targeting site that would impact the coupling between HL1/TM2 and the catalytic pocket, opening, to our knowledge, new avenues for structure-based design of novel allosteric modulators of γ-secretase protease activity.