Polymorph sampling with coupling to extended variables: enhanced sampling of polymorph energy landscapes and free energy perturbation of polymorph ensembles.

A novel approach to computationally enhance the sampling of molecular crystal structures is proposed and tested. This method is based on the use of extended variables coupled to a Monte Carlo based crystal polymorph generator. Inspired by the established technique of quasi-random sampling of polymorphs using the rigid molecule constraint, this approach represents molecular clusters as extended variables within a thermal reservoir. Polymorph unit-cell variables are generated using pseudo-random sampling. Within this framework, a harmonic coupling between the extended variables and polymorph configurations is established. The extended variables remain fixed during the inner loop dedicated to polymorph sampling, enforcing a stepwise propagation of the extended variables to maintain system exploration. The final processing step results in a polymorph energy landscape, where the raw structures sampled to create the extended variable trajectory are re-optimized without the thermal coupling term. The foundational principles of this approach are described and its effectiveness using both a Metropolis Monte Carlo type algorithm and modifications that incorporate replica exchange is demonstrated. A comparison is provided with pseudo-random sampling of polymorphs for the molecule coumarin. The choice to test a design of this algorithm as relevant for enhanced sampling of crystal structures was due to the obvious relation between molecular structure variables and corresponding crystal polymorphs as representative of the inherent vapor to crystal transitions that exist in nature. Additionally, it is shown that the trajectories of extended variables can be harnessed to extract fluctuation properties that can lead to valuable insights. A novel thermodynamic variable is introduced: the free energy difference between ensembles of Z' = 1 and Z' = 2 crystal polymorphs.

Adaptive CVgen: Leveraging reinforcement learning for advanced sampling in protein folding and chemical reactions.

Enhanced sampling techniques have traditionally encountered two significant challenges: identifying suitable reaction coordinates and addressing the exploration-exploitation dilemma, particularly the difficulty of escaping local energy minima. Here, we introduce Adaptive CVgen, a universal adaptive sampling framework designed to tackle these issues. Our approach utilizes a set of collective variables (CVs) to comprehensively cover the system's potential evolutionary phase space, generating diverse reaction coordinates to address the first challenge. Moreover, we integrate reinforcement learning strategies to dynamically adjust the generated reaction coordinates, thereby effectively balancing the exploration-exploitation dilemma. We apply this framework to sample the conformational space of six proteins transitioning from completely disordered states to folded states, as well as to model the chemical synthesis process of C60, achieving conformations that perfectly match the standard C60 structure. The results demonstrate Adaptive CVgen's effectiveness in exploring new conformations and escaping local minima, achieving both sampling efficiency and exploration accuracy. This framework holds potential for extending to various related challenges, including protein folding dynamics, drug targeting, and complex chemical reactions, thereby opening promising avenues for application in these fields.

Gating residues govern ligand unbinding kinetics from the buried cavity in HIF-2α PAS-B.

While transcription factors have been generally perceived as "undruggable," an exception is the HIF-2 hypoxia-inducible transcription factor, which contains an internal cavity that is sufficiently large to accommodate a range of small-molecules, including the therapeutically used inhibitor belzutifan. Given the relatively long ligand residence times of these small molecules and the lack of any experimentally observed pathway connecting the cavity to solvent, there has been great interest in understanding how these drug ligands exit the buried receptor cavity. Here, we focus on the relevant PAS-B domain of hypoxia-inducible factor 2α (HIF-2α) and examine how one such small molecule (THS-017) exits from the buried cavity within this domain on the seconds-timescale using atomistic simulations and ZZ-exchange NMR. To enable the simulations, we applied the weighted ensemble path sampling strategy, which generates continuous pathways for a rare-event process [e.g., ligand (un)binding] with rigorous kinetics in orders of magnitude less computing time compared to conventional simulations. Results reveal the formation of an encounter complex intermediate and two distinct classes of pathways for ligand exit. Based on these pathways, we identified two pairs of conformational gating residues in the receptor: one for the major class (N288 and S304) and another for the minor class (L272 and M309). ZZ-exchange NMR validated the kinetic importance of N288 for ligand unbinding. Our results provide an ideal simulation dataset for rational manipulation of ligand unbinding kinetics.

Empowering AlphaFold2 for protein conformation selective drug discovery with AlphaFold2-RAVE.

Small-molecule drug design hinges on obtaining co-crystallized ligand-protein structures. Despite AlphaFold2's strides in protein native structure prediction, its focus on apo structures overlooks ligands and associated holo structures. Moreover, designing selective drugs often benefits from the targeting of diverse metastable conformations. Therefore, direct application of AlphaFold2 models in virtual screening and drug discovery remains tentative. Here, we demonstrate an AlphaFold2-based framework combined with all-atom enhanced sampling molecular dynamics and Induced Fit docking, named AF2RAVE-Glide, to conduct computational model-based small-molecule binding of metastable protein kinase conformations, initiated from protein sequences. We demonstrate the AF2RAVE-Glide workflow on three different mammalian protein kinases and their type I and II inhibitors, with special emphasis on binding of known type II kinase inhibitors which target the metastable classical DFG-out state. These states are not easy to sample from AlphaFold2. Here, we demonstrate how with AF2RAVE these metastable conformations can be sampled for different kinases with high enough accuracy to enable subsequent docking of known type II kinase inhibitors with more than 50% success rates across docking calculations. We believe the protocol should be deployable for other kinases and more proteins generally.

Mechanisms of peptide agonist dissociation and deactivation of adhesion G-protein-coupled receptors.

Adhesion G protein-coupled receptors (ADGRs) belong to Class B2 of GPCRs and are involved in a wide array of important physiological processes. ADGRs contain a GPCR autoproteolysis-inducing (GAIN) domain that is proximal to the receptor N-terminus and undergoes autoproteolysis during biosynthesis to generate two fragments: the N-terminal fragment (NTF) and C-terminal fragment (CTF). Dissociation of NTF reveals a tethered agonist to activate CTF of ADGRs for G protein signaling. Synthetic peptides that mimic the tethered agonist can also activate the ADGRs. However, mechanisms of peptide agonist dissociation and deactivation of ADGRs remain poorly understood. In this study, we have performed all-atom enhanced sampling simulations using a novel Protein-Protein Interaction-Gaussian accelerated Molecular Dynamics (PPI-GaMD) method on the ADGRG2-IP15 and ADGRG1-P7 complexes. The PPI-GaMD simulations captured dissociation of the IP15 and P7 peptide agonists from their target receptors. We were able to identify important low-energy conformations of ADGRG2 and ADGRG1 in the active, intermediate, and inactive states, as well as exploring different states of the peptide agonists IP15 and P7 during dissociation. Therefore, our PPI-GaMD simulations have revealed dynamic mechanisms of peptide agonist dissociation and deactivation of ADGRG1 and ADGRG2, which will facilitate rational design of peptide regulators of the two receptors and other ADGRs.

Hidden Structural States of Proteins Revealed by Conformer Selection with AlphaFold-NMR.

Recent advances in molecular modeling using deep learning can revolutionize our understanding of dynamic protein structures. NMR is particularly well-suited for determining dynamic features of biomolecular structures. The conventional process for determining biomolecular structures from experimental NMR data involves its representation as conformation-dependent restraints, followed by generation of structural models guided by these spatial restraints. Here we describe an alternative approach: generating a distribution of realistic protein conformational models using artificial intelligence-(AI-) based methods and then selecting the sets of conformers that best explain the experimental data. We applied this conformational selection approach to redetermine the solution NMR structure of the enzyme Gaussia luciferase. First, we generated a diverse set of conformer models using AlphaFold2 (AF2) with an enhanced sampling protocol. The models that best-fit NOESY and chemical shift data were then selected with a Bayesian scoring metric. The resulting models include features of both the published NMR structure and the standard AF2 model generated without enhanced sampling. This "AlphaFold-NMR" protocol also generated an alternative "open" conformational state that fits nearly as well to the overall NMR data but accounts for some NOESY data that is not consistent with first "closed" conformational state; while other NOESY data consistent with this second state are not consistent with the first conformational state. The structure of this "open" structural state differs from that of the "closed" state primarily by the position of a thumb-shaped loop between α-helices H5 and H6, revealing a cryptic surface pocket. These alternative conformational states of Gluc are supported by "double recall" analysis of NOESY data and AF2 models. Additional structural states are also indicated by backbone chemical shift data indicating partially-disordered conformations for the C-terminal segment. Considered as a multistate ensemble, these multiple states of Gluc together fit the NOESY and chemical shift data better than the "restraint-based" NMR structure and provide novel insights into its structure-dynamic-function relationships. This study demonstrates the potential of AI-based modeling with enhanced sampling to generate conformational ensembles followed by conformer selection with experimental data as an alternative to conventional restraint satisfaction protocols for protein NMR structure determination.

Molecular Dynamics Simulations in Protein-Protein Docking.

Concerted interactions between all the cell components form the basis of biological processes. Protein-protein interactions (PPIs) constitute a tremendous part of this interaction network. Deeper insight into PPIs can help us better understand numerous diseases and lead to the development of new diagnostic and therapeutic strategies. PPI interfaces, until recently, were considered undruggable. However, it is now believed that the interfaces contain "hot spots," which could be targeted by small molecules. Such a strategy would require high-quality structural data of PPIs, which are difficult to obtain experimentally. Therefore, in silico modeling can complement or be an alternative to in vitro approaches. There are several computational methods for analyzing the structural data of the binding partners and modeling of the protein-protein dimer/oligomer structure. The major problem with in silico structure prediction of protein assemblies is obtaining sufficient sampling of protein dynamics. One of the methods that can take protein flexibility and the effects of the environment into account is Molecular Dynamics (MD). While sampling of the whole protein-protein association process with plain MD would be computationally expensive, there are several strategies to harness the method to PPI studies while maintaining reasonable use of resources. This chapter reviews known applications of MD in the PPI investigation workflows.

Design, synthesis, pharmacological evaluation, and computational study of benzo[d] isothiazol-based small molecule inhibitors targeting PD-1/PD-l1 interaction.

Blockade of the programmed cell death-1 (PD-1)/programmed cell death ligand 1 (PD-L1) pathway is an attractive strategy for immunotherapy, but the clinical application of small molecule PD-1/PD-L1 inhibitors remains unclear. In this work, based on BMS-202 and our previous work YLW-106, a series of compounds with benzo[d]isothiazol structure as scaffold were designed and synthesized. Their inhibitory activity against PD-1/PD-L1 interaction was evaluated by a homogeneous time-resolved fluorescence (HTRF) assay. Among them, LLW-018 (27c) exhibited the most potent inhibitory activity with an IC50 value of 2.61 nM. The cellular level assays demonstrated that LLW-018 exhibited low cytotoxicity against Jurkat T and MDA-MB-231. Further cell-based PD-1/PD-L1 blockade bioassays based on PD-1 NFAT-Luc Jurkat cells and PD-L1 TCR Activator CHO cells indicated that LLW-018 could interrupt PD-1/PD-L1 interaction with an IC50 value of 0.88 µM. Multi-computational methods, including molecular docking, molecular dynamics, MM/GBSA, MM/PBSA, Metadynamics, and QM/MM MD were utilized on PD-L1 dimer complexes, which revealed the binding modes and dissociation process of LLW-018 and C2-symmetric small molecule inhibitor LCH1307. These results suggested that LLW-018 exhibited promising potency as a PD-1/PD-L1 inhibitor for further investigation.

Dynamical Reweighting for Biased Rare Event Simulations.

Dynamical reweighting techniques aim to recover the correct molecular dynamics from a simulation at a modified potential energy surface. They are important for unbiasing enhanced sampling simulations of molecular rare events. Here, we review the theoretical frameworks of dynamical reweighting for modified potentials. Based on an overview of kinetic models with increasing level of detail, we discuss techniques to reweight two-state dynamics, multistate dynamics, and path integrals. We explore the natural link to transition path sampling and how the effect of nonequilibrium forces can be reweighted. We end by providing an outlook on how dynamical reweighting integrates with techniques for optimizing collective variables and with modern potential energy surfaces.

Molecular Mechanisms of the Impaired Heparin Pentasaccharide Interactions in 10 Antithrombin Heparin Binding Site Mutants Revealed by Enhanced Sampling Molecular Dynamics.

Antithrombin (AT) is a critical regulator of the coagulation cascade by inhibiting multiple coagulation factors including thrombin and FXa. Binding of heparinoids to this serpin enhances the inhibition considerably. Mutations located in the heparin binding site of AT result in thrombophilia in affected individuals. Our aim was to study 10 antithrombin mutations known to affect their heparin binding in a heparin pentasaccharide bound state using two molecular dynamics (MD) based methods providing enhanced sampling, GaMD and LiGaMD2. The latter provides an additional boost to the ligand and the most important binding site residues. From our GaMD simulations we were able to identify four variants (three affecting amino acid Arg47 and one affecting Lys114) that have a particularly large effect on binding. The additional acceleration provided by LiGaMD2 allowed us to study the consequences of several other mutants including those affecting Arg13 and Arg129. We were able to identify several conformational types by cluster analysis. Analysis of the simulation trajectories revealed the causes of the impaired pentasaccharide binding including pentasaccharide subunit conformational changes and altered allosteric pathways in the AT protein. Our results provide insights into the effects of AT mutations interfering with heparin binding at an atomic level and can facilitate the design or interpretation of in vitro experiments.

Application of Funnel Metadynamics to the Platelet Integrin αIIbß3 in Complex with an RGD Peptide.

Integrin αIIbß3 mediates platelet aggregation by binding the Arginyl-Glycyl-Aspartic acid (RGD) sequence of fibrinogen. RGD binding occurs at a site topographically proximal to the αIIb and ß3 subunits, promoting the conformational activation of the receptor from bent to extended states. While several experimental approaches have characterized RGD binding to αIIbß3 integrin, applying computational methods has been significantly more challenging due to limited sampling and the need for a priori information regarding the interactions between the RGD peptide and integrin. In this study, we employed all-atom simulations using funnel metadynamics (FM) to evaluate the interactions of an RGD peptide with the αIIb and ß3 subunits of integrin. FM incorporates an external history-dependent potential on selected degrees of freedom while applying a funnel-shaped restraint potential to limit RGD exploration of the unbound state. Furthermore, it does not require a priori information about the interactions, enhancing the sampling at a low computational cost. Our FM simulations reveal significant molecular changes in the ß3 subunit of integrin upon RGD binding and provide a free-energy landscape with a low-energy binding mode surrounded by higher-energy prebinding states. The strong agreement between previous experimental and computational data and our results highlights the reliability of FM as a method for studying dynamic interactions of complex systems such as integrin.

Ligand Gaussian accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides.

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional Molecular Dynamics (cMD), due to limited simulation timescales. Based on our previously developed Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3", in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 microsecond simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 agreed with available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligand and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

Unlocking the specificity of antimicrobial peptide interactions for membrane-targeted therapies.

Antimicrobial peptides (AMPs) are increasingly recognized as potent therapeutic agents, with their selective affinity for pathological membranes, low toxicity profile, and minimal resistance development making them particularly attractive in the pharmaceutical landscape. This study offers a comprehensive analysis of the interaction between specific AMPs, including magainin-2, pleurocidin, CM15, LL37, and clavanin, with lipid bilayer models of very different compositions that have been ordinarily used as biological membrane models of healthy mammal, cancerous, and bacterial cells. Employing unbiased molecular dynamics simulations and metadynamics techniques, we have deciphered the intricate mechanisms by which these peptides recognize pathogenic and pathologic lipid patterns and integrate into lipid assemblies. Our findings reveal that the transverse component of the peptide's hydrophobic dipole moment is critical for membrane interaction, decisively influencing the molecule's orientation and expected therapeutic efficacy. Our approach also provides insight on the kinetic and dynamic dependence on the peptide orientation in the axial and azimuthal angles when coming close to the membrane. The aim is to establish a robust framework for the rational design of peptide-based, membrane-targeted therapies, as well as effective quantitative descriptors that can facilitate the automated design of novel AMPs for these therapies using machine learning methods.

Understanding the impact of binding free energy and kinetics calculations in modern drug discovery.

INTRODUCTION: For rational drug design, it is crucial to understand the receptor-drug binding processes and mechanisms. A new era for the use of computer simulations in predicting drug-receptor interactions at an atomic level has begun with remarkable advances in supercomputing and methodological breakthroughs. AREAS COVERED: End-point free energy calculation methods such as Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) or Molecular-Mechanics/Generalized Born Surface Area (MM/GBSA), free energy perturbation (FEP), and thermodynamic integration (TI) are commonly used for binding free energy calculations in drug discovery. In addition, kinetic dissociation and association rate constants (koff and kon) play critical roles in the function of drugs. Nowadays, Molecular Dynamics (MD) and enhanced sampling simulations are increasingly being used in drug discovery. Here, the authors provide a review of the computational techniques used in drug binding free energy and kinetics calculations. EXPERT OPINION: The applications of computational methods in drug discovery and design are expanding, thanks to improved predictions of the binding free energy and kinetic rates of drug molecules. Recent microsecond-timescale enhanced sampling simulations have made it possible to accurately capture repetitive ligand binding and dissociation, facilitating more efficient and accurate calculations of ligand binding free energy and kinetics.

[Development of Membrane Permeability Coefficient by Means of Novel Molecular Dynamics Methods].

The membrane permeability, and its evaluation, is crucial factor in the process of uptake of compounds from outside to inside the cell and in the inhibition of the activity of disease-causing target proteins. Although molecular dynamics (MD) simulations have been shown to be able to reproduce the conformational changes of compounds occurring during membrane permeation, it is still challenging to extract the membrane permeability at an affordable computational workload solely by conventional MD. Indeed, the time scale accessible by MD is far below the one characterizing the actual permeation process. Phenomena occurring in living organisms escaping the reach of standard MD are generally referred to as biological rare events, and the membrane permeation process is one of them. To overcome this time-scale problem, several enhanced sampling methods have been proposed over the years to improve conformational sampling. In this review, a hybrid sampling method that combines the parallel cascade selection MD (PaCS-MD) and the outlier flooding method (OFLOOD), introduced and developed by our group, is proposed as a tool to study the membrane permeation from structural sampling (rare-event sampling). The obtained trajectories are used to estimate the free energy profiles for the membrane permeation and to compute the membrane permeation coefficients. Moreover, we present an example of application of the free energy reaction network method as a versatile way for incorporating explicitly into reaction coordinates the degrees of freedom related to internal motion.

Effective data-driven collective variables for free energy calculations from metadynamics of paths.

A variety of enhanced sampling (ES) methods predict multidimensional free energy landscapes associated with biological and other molecular processes as a function of a few selected collective variables (CVs). The accuracy of these methods is crucially dependent on the ability of the chosen CVs to capture the relevant slow degrees of freedom of the system. For complex processes, finding such CVs is the real challenge. Machine learning (ML) CVs offer, in principle, a solution to handle this problem. However, these methods rely on the availability of high-quality datasets-ideally incorporating information about physical pathways and transition states-which are difficult to access, therefore greatly limiting their domain of application. Here, we demonstrate how these datasets can be generated by means of ES simulations in trajectory space via the metadynamics of paths algorithm. The approach is expected to provide a general and efficient way to generate efficient ML-based CVs for the fast prediction of free energy landscapes in ES simulations. We demonstrate our approach with two numerical examples, a 2D model potential and the isomerization of alanine dipeptide, using deep targeted discriminant analysis as our ML-based CV of choice.

The Recognition Pathway of the SARS-CoV-2 Spike Receptor-Binding Domain to Human Angiotensin-Converting Enzyme 2.

COVID-19 caused by SARS-CoV-2 has spread around the world. The receptor-binding domain (RBD) of the spike protein of SARS-CoV-2 is a critical component that directly interacts with host ACE2. Here, we simulate the ACE2 recognition processes of RBD of the WT, Delta, and OmicronBA.2 variants using our recently developed supervised Gaussian accelerated molecular dynamics (Su-GaMD) approach. We show that RBD recognizes ACE2 through three contact regions (regions I, II, and III), which aligns well with the anchor-locker mechanism. The higher binding free energy in State d of the RBDOmicronBA.2-ACE2 system correlates well with the increased infectivity of OmicronBA.2 in comparison with other variants. For RBDDelta, the T478K mutation affects the first step of recognition, while the L452R mutation, through its nearby Y449, affects the RBDDelta-ACE2 binding in the last step of recognition. For RBDOmicronBA.2, the E484A mutation affects the first step of recognition, the Q493R, N501Y, and Y505H mutations affect the binding free energy in the last step of recognition, mutations in the contact regions affect the recognition directly, and other mutations indirectly affect recognition through dynamic correlations with the contact regions. These results provide theoretical insights for RBD-ACE2 recognition and may facilitate drug design against SARS-CoV-2.

Insights into the mechanism of peptide fibril growth on gold surface.

Understanding the formation of ß-fibrils over the gold surface is of paramount interest in nano-bio-medicinal Chemistry. The intricate mechanism of self-assembly of neurofibrillogenic peptides and their growth over the gold surface remains elusive, as experiments are limited in unveiling the microscopic dynamic details, in particular, at the early stage of the peptide aggregation. In this work, we carried out equilibrium molecular dynamics and enhanced sampling simulations to elucidate the underlying mechanism of the growth of an amyloid-forming sequence of tau fragments over the gold surface. Our results disclose that the collective intermolecular interactions between the peptide chains and peptides with the gold surface facilitate the peptide adsorption, followed by integration, finally leading to the fibril formation.

Investigation of the MDM2-binding potential of de novo designed peptides using enhanced sampling simulations.

The tumor suppressor p53 plays a crucial role in cellular responses to various stresses, regulating key processes such as apoptosis, senescence, and DNA repair. Dysfunctional p53, prevalent in approximately 50 % of human cancers, contributes to tumor development and resistance to treatment. This study employed deep learning-based protein design and structure prediction methods to identify novel high-affinity peptide binders (Pep1 and Pep2) targeting MDM2, with the aim of disrupting its interaction with p53. Extensive all-atom molecular dynamics simulations highlighted the stability of the designed peptide in complex with the target, supported by several structural analyses, including RMSD, RMSF, Rg, SASA, PCA, and free energy landscapes. Using the steered molecular dynamics and umbrella sampling simulations, we elucidate the dissociation dynamics of p53, Pep1, and Pep2 from MDM2. Notable differences in interaction profiles were observed, emphasizing the distinct dissociation patterns of each peptide. In conclusion, the results of our umbrella sampling simulations suggest Pep1 as a higher-affinity MDM2 binder compared to p53 and Pep2, positioning it as a potential inhibitor of the MDM2-p53 interaction. Using state-of-the-art protein design tools and advanced MD simulations, this study provides a comprehensive framework for rational in silico design of peptide binders with therapeutic implications in disrupting MDM2-p53 interactions for anticancer interventions.

Enhanced Sampling with Machine Learning.

Molecular dynamics (MD) enables the study of physical systems with excellent spatiotemporal resolution but suffers from severe timescale limitations. To address this, enhanced sampling methods have been developed to improve the exploration of configurational space. However, implementing these methods is challenging and requires domain expertise. In recent years, integration of machine learning (ML) techniques into different domains has shown promise, prompting their adoption in enhanced sampling as well. Although ML is often employed in various fields primarily due to its data-driven nature, its integration with enhanced sampling is more natural with many common underlying synergies. This review explores the merging of ML and enhanced MD by presenting different shared viewpoints. It offers a comprehensive overview of this rapidly evolving field, which can be difficult to stay updated on. We highlight successful strategies such as dimensionality reduction, reinforcement learning, and flow-based methods. Finally, we discuss open problems at the exciting ML-enhanced MD interface.