Enhanced Sampling Molecular Dynamics Research Articles

Stabilized SARS-CoV-2 spike S2 immunogens designed using insights from enhanced sampling molecular dynamics simulations.

The SARS-CoV-2 spike protein is an important molecular machine that facilitates membrane fusion to host cells and is the major antigenic component of vaccines. Notably, full-length prefusion-stabilized spike proteins produce a strong neutralizing antibody response against S1, which contains the N-terminal domain (NTD) and receptor binding domain (RBD). Unfortunately, emerging variants of concern contain numerous mutations, insertions, and deletions within S1 that impact the binding of neutralizing antibodies which may diminish the efficacy of current vaccines. In contrast, S2 is the most conserved region of the spike and is capable of eliciting antibodies that bind broadly to diverse β-coronavirus spikes. Although S2-based immunogens may have potential benefits, a challenge to its usage as an alternative vaccine strategy is the general instability of S2 in the absence of S1. Using enhanced sampling molecular dynamics simulations, we observed that the S2 head adopts an open conformation where individual protomers separate from one another at the apex. Extensive contact map analysis reveals that the formation and breakage of contacts between three helices mediates these opening motions. Through introducing mutations in these helices, we designed a stabilized S2 construct that maintains a closed conformation and has enhanced thermal stability and protein expression. These studies represent the first simulations characterizing S2 conformational dynamics and demonstrate the potential of using enhanced sampling molecular dynamics simulations in the design of stabilized SARS-CoV-2 S2 immunogens.

PIP2 binding at the voltage sensor domain facilitates KCNQ1 VSD activation and gating.

The cardiac KCNQ1 channel is involved in cell hyperpolarisation, and loss-of-function mutations are associated with Long QT syndrome (LQTS), which can cause sudden cardiac arrest. The KCNQ1 channel is gated by the changes in transmembrane voltage and crucially requires the binding of two phosphatidylinositol-4,5-bisphosphate (PIP2) molecules to open the channel. The recent structure of KCNQ1 has revealed a PIP2 binding site at the voltage-sensor domain (VSD). In addition, previous modelling and electrophysiological studies have proposed another PIP2 binding site at the pore domain. Using enhanced sampling molecular dynamics simulation, we show that PIP2 binding at the VSD facilitates its activation, whereas PIP2 interaction of the pore facilitates VSD-pore coupling by straightening the S6 helix. Our simulation proposes that PIP2 binding at the VSD acts as a counter charge pausing point during the voltage-dependent activation process, lowering the energy barrier during the S4 translocation. Together, our work highlights the importance of PIP2 at the VSD as an activator to the VSD activation, in addition to our previous knowledge as the key coupler between the VSD and the pore.

Mutations in a human kinase cause drug resistance by shifting conformational state distributions.

In pharmacology, protein structures are frequently represented as static snapshots with mostly rigid conformations. This portrayal is likely myopic, as both theory and experiment have shown that biological macromolecules tend to occupy time-dependent conformational ensembles with thousands of states, each of which can adopt distinct functional and binding profiles. We believe that studying the often disregarded dynamics of these macromolecular structures should take center stage in the next generation of drug discovery efforts. This philosophy stems from the intuition that ensembles and their internal transitions multiply the number of potentially druggable targets. Recognizing this paradigm shift, we have deployed enhanced-sampling molecular dynamics (MD) simulations to map the conformational landscape of cancer-associated protein kinase domains. Through MD, we have successfully sampled the transition from the active state (known as DFG in) of wild-type and mutant forms of the Abelson tyrosine kinase (Abl1) to two distinct inactive states: DFG out and DFG out/collapsed A-loop. This transition spans an overall free energy barrier of over 80 kcal/mol from the active state to the second inactive state. Our work has revealed a rugged conformational landscape, with a collection of potentially metastable conformations that are virtually undetectable with traditional methods and could be used as targets for drug discovery. Importantly, we have also derived molecular mechanisms for the observed transitions between states, with key implications for biochemistry and protein engineering. Finally, we observed that two mutations that are known to cause inhibitor resistance in cancer patients lead to their phenotypes by significantly shifting the free energy barriers and pathways between states, locking the kinases in their active conformations. Ultimately, our work could lead to exciting opportunities for pharmacology, and further cement the importance of macromolecular conformational dynamics in biomedical discovery.

Molecular determinants of distinctive opioid receptor subtype affinities.

Despite being targeted by over 30% of the drugs available on the market, G-Protein Coupled Receptors (GPCR) have proven to be notoriously difficult to characterize on an atomistic level. While it is known that helix 6 (TM6) movement is necessary for signal transduction, structural details of activation remain elusive and each receptor is generally studied individually, making comparisons between homologous proteins challenging. Given their homology, any newly discovered drug targeting a specific GPCR may be promiscuous, adding work- and cost-intensive dimensions to the drug discovery process. Some of the GPCRs that seem the most appealing as drug targets are opioid receptors. As opioids are potent pain relievers with widespread use and addictive properties, understanding how drugs bind particular opioid receptor subtypes would allow the design of non-addictive derivatives and thus help alleviate the ongoing opioid crisis - in the United States only, it is estimated that 4 million people use opioids recreationally or are dependent on them, with about half overdose-related deaths being related to synthetic opioids. In this work, we have used enhanced-sampling molecular dynamics simulations to explore the conformational dynamics of different subtypes of opioid receptors and uncovered the molecular determinants of their affinity to various endogenous peptides and synthetic opioids. Using enhanced sampling techniques, we have obtained free energy surfaces (FES) of activation of the μ-, κ-, and δ-subtypes. By comparing the metastable states obtained from the FES, we have pinpointed key structural differences relevant to the activation process. These findings will allow the development of simple and accurate structure-activity relationship, which in turn will guide the design of new, less addictive opioids.

Using tetrapeptide GSSS and β-catenin17-48 to study the effects of phosphoserine on polypeptide conformational ensembles.

Phosphorylation is one of the most common post-translational modifications and is known to have roles in protein transport, localization, function, and stability. Studying both the phosphorylated and nonphosphorylated states of a protein is important as they can be correlated to health and disease states, such as in Alzheimer's and Parkinson's disease. β-catenin is a protein that has two independent functions: cell-cell adhesion and transcriptional regulation, and when there is a change in the balance of these roles, β-catenin acts as an oncogene. It is thought that the phosphorylation state of the unstructured N-terminal region plays an important role in regulating the function and localization of β-catenin. Additionally, mutations at Ser33 and Ser37 phosphorylation sites have been associated with several types of cancer including liver, ovarian, and colorectal cancer. The addition of a large, negatively charged group on the serine sidechain is expected to strongly perturb local electric fields, thus influencing backbone conformational sampling and propensity for hydrogen bonding. As such, dipole-dipole interactions are expected to feature prominently in these systems and are best suited to be investigated with an explicitly polarizable force field. Parameters for phosphoserine were recently developed for the Drude polarizable force field. To validate these parameters and study a biologically relevant system, we performed both unbiased and enhanced-sampling molecular dynamics simulations of the phosphorylated and nonphosphorylated states of model tetrapeptide GSSS and a peptide derived from the N-terminal region of β-catenin, including residues 17-48 (β-catenin17-48). We calculated 3JHNHα-couplings to compare against experimental data and assessed the effect of phosphorylation on secondary structure content, peptide-bond dipole moments, and the electric fields acting around the phosphorylation sites. In doing so, we seek to relate electronic structure with polypeptide conformational dynamics.

MDexciteR: Enhanced Sampling Molecular Dynamics by Excited Normal Modes or Principal Components Obtained from Experiments.

Molecular dynamics with excited normal modes (MDeNM) is an enhanced sampling method for exploring conformational changes in proteins with minimal biases. The excitation corresponds to injecting kinetic energy along normal modes describing intrinsic collective motions. Herein, we developed a new automated open-source implementation, MDexciteR (https://github.com/mcosta27/MDexciteR), enabling the integration of MDeNM with two commonly used simulation programs with GPU support. Second, we generalized the method to include the excitation of principal components calculated from experimental ensembles. Finally, we evaluated whether the use of coarse-grained normal modes calculated with elastic network representations preserved the performance and accuracy of the method. The advantages and limitations of these new approaches are discussed based on results obtained for three different protein test cases: two globular and a protein/membrane system.

Dynamical characterization and multiple unbinding paths of two PreQ1 ligands in one pocket.

Riboswitches naturally regulate gene expression in bacteria by binding to specific small molecules. Class 1 preQ1 riboswitch aptamer is an important model not only for RNA folding but also as a target for designing small molecule antibiotics due to its well-known minimal aptamer domain. Here, we ran a total of 62.4 μs conventional and enhanced-sampling molecular dynamics (MD) simulations to characterize the determinants underlying the binding of the preQ1-II riboswitch aptamer to two preQ1 ligands in one binding pocket. Decomposition of binding free energy suggested that preQ1 ligands at α and β sites interact with four nucleotides (G5, C17, C18, and A30) and two nucleotides (A12 and C31), respectively. Mg2+ ions play a crucial role in both stabilizing the binding pocket and facilitating ligand binding. The flexible preQ1 ligand at the β site leads to the top of the binding pocket loosening and thus pre-organizes the riboswitch for ligand entry. Enhanced sampling simulations further revealed that the preQ1 ligand at the α site unbinds through two orthogonal pathways, which are dependent on whether or not a β site preQ1 ligand is present. One of the two preQ1 ligands has been identified in the binding pocket, which will aid to identify the second preQ1 Ligand. Our work provides new information for designing robust ligands.

Tidying up the conformational ensemble of a disordered peptide by computational prediction of spectroscopic fingerprints.

The most advanced structure prediction methods are powerless in exploring the conformational ensemble of disordered peptides and proteins and for this reason the "protein folding problem" remains unsolved. We present a novel methodology that enables the accurate prediction of spectroscopic fingerprints (circular dichroism, infrared, Raman, and Raman optical activity), and by this allows for "tidying up" the conformational ensembles of disordered peptides and disordered regions in proteins. This concept is elaborated for and applied to a dodecapeptide, whose spectroscopic fingerprint is measured and theoretically predicted by means of enhanced-sampling molecular dynamics coupled with quantum mechanical calculations. Following this approach, we demonstrate that peptides lacking a clear propensity for ordered secondary-structure motifs are not randomly, but only conditionally disordered. This means that their conformational landscape, or phase-space, can be well represented by a basis-set of conformers including about 10 to 100 structures. The implications of this finding have profound consequences both for the interpretation of experimental electronic and vibrational spectral features of peptides in solution and for the theoretical prediction of these features using accurate and computationally expensive techniques. The here-derived methods and conclusions are expected to fundamentally impact the rationalization of so-far elusive structure-spectra relationships for disordered peptides and proteins, towards improved and versatile structure prediction methods.

A combined ligand and target-based virtual screening strategy to repurpose drugs as putrescine uptake inhibitors with trypanocidal activity.

Chagas disease, also known as American trypanosomiasis, is a neglected tropical disease caused by the protozoa Trypanosoma cruzi, affecting nearly 7million people only in the Americas. Polyamines are essential compounds for parasite growth, survival, and differentiation. However, because trypanosomatids are auxotrophic for polyamines, they must be obtained from the host by specific transporters. In this investigation, an ensemble of QSAR classifiers able to identify polyamine analogs with trypanocidal activity was developed. Then, a multi-template homology model of the dimeric polyamine transporter of T. cruzi, TcPAT12, was created with Rosetta, and then refined by enhanced sampling molecular dynamics simulations. Using representative snapshots extracted from the trajectory, a docking model able to discriminate between active and inactive compounds was developed and validated. Both models were applied in a parallel virtual screening campaign to repurpose known drugs as anti-trypanosomal compounds inhibiting polyamine transport in T. cruzi. Montelukast, Quinestrol, Danazol, and Dutasteride were selected for in vitro testing, and all of them inhibited putrescine uptake in biochemical assays, confirming the predictive ability of the computational models. Furthermore, all the confirmed hits proved to inhibit epimastigote proliferation, and Quinestrol and Danazol were able to inhibit, in the low micromolar range, the viability of trypomastigotes and the intracellular growth of amastigotes.

Kinetics of Drug Release from Clay Using Enhanced Sampling Methods.

A key step in the development of a new drug, is the design of drug-excipient complexes that lead to optimal drug release kinetics. Computational chemistry and specifically enhanced sampling molecular dynamics methods can play a key role in this context, by minimizing the need for expensive experiments, and reducing cost and time. Here we show that recent advances in enhanced sampling methodologies can be brought to fruition in this area. We demonstrate the potential of these methodologies by simulating the drug release kinetics of the complex praziquantel-montmorillonite in water. Praziquantel finds promising applications in the treatment of schistosomiasis, but its biopharmaceutical profile needs to be improved, and a cheap material such as the montmorillonite clay would be a very convenient excipient. We simulate the drug release both from surface and interlayer space, and find that the diffusion of the praziquantel inside the interlayer space is the process that limits the rate of drug release.

Non-uniform Self-folding of Helical Poly(fluorenevinylene) Derivatives in the Solid State Leading to Amplified Circular Dichroism and Circularly Polarized Light Emission.

An unprecedented non-uniform self-folding of artificial polymer chains composed of turn moieties and stretched segments is presented through the design of a set of optically active poly(fluorene-2,7-diylethene-1,2-diyl) (poly(fluorenevinylene)) derivatives bearing a neomenthyl group and a pentyl group attached at the 9-position of fluorene backbone at various ratios. The folded structure is formed and stabilized through inter-chain interactions in the solid state, leading to remarkably enhanced chiroptical properties (chirality amplification) in terms of circular dichroism (CD) and circularly polarized light (CPL) emission. This phenomenon is rationalized by experimental and theoretical CD and CPL spectral analyses. The polymer arrangements in the solid state were further assessed through transmission electron microscopic observations combined with enhanced sampling molecular dynamics simulations in the solid state revealing the thin film organizations.

Non‐uniform Self‐folding of Helical Poly(fluorenevinylene) Derivatives in the Solid State Leading to Amplified Circular Dichroism and Circularly Polarized Light Emission

AbstractAn unprecedented non‐uniform self‐folding of artificial polymer chains composed of turn moieties and stretched segments is presented through the design of a set of optically active poly(fluorene‐2,7‐diylethene‐1,2‐diyl) (poly(fluorenevinylene)) derivatives bearing a neomenthyl group and a pentyl group attached at the 9‐position of fluorene backbone at various ratios. The folded structure is formed and stabilized through inter‐chain interactions in the solid state, leading to remarkably enhanced chiroptical properties (chirality amplification) in terms of circular dichroism (CD) and circularly polarized light (CPL) emission. This phenomenon is rationalized by experimental and theoretical CD and CPL spectral analyses. The polymer arrangements in the solid state were further assessed through transmission electron microscopic observations combined with enhanced sampling molecular dynamics simulations in the solid state revealing the thin film organizations.

Modified Protein-Water Interactions in CHARMM36m for Thermodynamics and Kinetics of Proteins in Dilute and Crowded Solutions.

Proper balance between protein-protein and protein-water interactions is vital for atomistic molecular dynamics (MD) simulations of globular proteins as well as intrinsically disordered proteins (IDPs). The overestimation of protein-protein interactions tends to make IDPs more compact than those in experiments. Likewise, multiple proteins in crowded solutions are aggregated with each other too strongly. To optimize the balance, Lennard-Jones (LJ) interactions between protein and water are often increased about 10% (with a scaling parameter, λ = 1.1) from the existing force fields. Here, we explore the optimal scaling parameter of protein-water LJ interactions for CHARMM36m in conjunction with the modified TIP3P water model, by performing enhanced sampling MD simulations of several peptides in dilute solutions and conventional MD simulations of globular proteins in dilute and crowded solutions. In our simulations, 10% increase of protein-water LJ interaction for the CHARMM36m cannot maintain stability of a small helical peptide, (AAQAA)3 in a dilute solution and only a small modification of protein-water LJ interaction up to the 3% increase (λ = 1.03) is allowed. The modified protein-water interactions are applicable to other peptides and globular proteins in dilute solutions without changing thermodynamic properties from the original CHARMM36m. However, it has a great impact on the diffusive properties of proteins in crowded solutions, avoiding the formation of too sticky protein-protein interactions.

Toward Computing Accurate Free Energies in Heterogeneous Catalysis: a Case Study for Adsorbed Isobutene in H-ZSM-5.

Herein, we propose a novel computational protocol that enables calculating free energies with improved accuracy by combining the best available techniques for enthalpy and entropy calculation. While the entropy is described by enhanced sampling molecular dynamics techniques, the energy is calculated using ab initio methods. We apply the method to assess the stability of isobutene adsorption intermediates in the zeolite H-SSZ-13, a prototypical problem that is computationally extremely challenging in terms of calculating enthalpy and entropy. We find that at typical operating conditions for zeolite catalysis (400 °C), the physisorbed π-complex, and not the tertiary carbenium ion as often reported, is the most stable intermediate. This method paves the way for sampling-based techniques to calculate the accurate free energies in a broad range of chemistry-related disciplines, thus presenting a big step forward toward predictive modeling.

Improving All-Atom Force Field to Accurately Describe DNA G-Quadruplex Loops.

The DNA G-quadruplex (GQ) displays structural polymorphisms, and interactions between its loops and flanking sequences critically determine which of the diverse GQ conformers is adopted. All-atom molecular dynamics (MD) simulations of GQs are computationally challenging due to slow folding times and force field (ff) artifacts. In an earlier study, a direct folding simulation of the simplest DNA GQ (TBA15) was first reported using a modified version of the AMBER bsc1 ff (bsc1_vdW ff). Despite this successful folding simulation, it was later found that the bsc1_vdW ff is somewhat limited in terms of describing loop structures of GQs, which is problematic because GQ loop regions play key roles in ligand binding to modulate GQ activities. In this study, we further modified the bsc1_vdW ff to enhance the GQ loop prediction by fine-tuning a limited number of van der Waals (vdW) parameters of the standard AMBER bsc1 ff to improve the GQ loop distribution of a target GQ system (three-layered antiparallel GQ; mHtel21). Test simulations of this newly generated ff (bsc1_vdWL ff) on DNA GQs with diverse topologies (hybrid1, hybrid2, and parallel propeller) revealed that loop structures were predicted more accurately than by the bsc1_vdW ff. We consider that enhanced sampling MD simulation methods in combination with bsc1_vdWL provide useful simulation protocols for resolving outstanding issues of DNA GQ folding and GQ/ligand binding at the all-atom level.

Phase diagram of the TIP4P/Ice water model by enhanced sampling simulations.

We studied the phase diagram for the TIP4P/Ice water model using enhanced sampling molecular dynamics simulations. Our approach is based on the calculation of ice-liquid free energy differences from biased coexistence simulations that reversibly sample the melting and growth of layers of ice. We computed a total of 19 melting points for five different ice polymorphs, which are in excellent agreement with the melting lines obtained from the integration of the Clausius-Clapeyron equation. For proton-ordered and fully proton-disordered ice phases, the results are in very good agreement with previous calculations based on thermodynamic integration. For the partially proton-disordered ice III, we find a large increase in stability that is in line with previous observations using direct coexistence simulations for the TIP4P/2005 model. This issue highlights the robustness of the approach employed here for ice polymorphs with diverse degrees of proton disorder. Our approach is general and can be applied to the calculation of other complex phase diagrams.

Conformational Entropy as a Potential Liability of Computationally Designed Antibodies.

In silico antibody discovery is emerging as a viable alternative to traditional in vivo and in vitro approaches. Many challenges, however, remain open to enabling the properties of designed antibodies to match those produced by the immune system. A major question concerns the structural features of computer-designed complementarity determining regions (CDRs), including the role of conformational entropy in determining the stability and binding affinity of the designed antibodies. To address this problem, we used enhanced-sampling molecular dynamics simulations to compare the free energy landscapes of single-domain antibodies (sdAbs) designed using structure-based (DesAb-HSA-D3) and sequence-based approaches (DesAbO), with that of a nanobody derived from llama immunization (Nb10). Our results indicate that the CDR3 of DesAbO is more conformationally heterogeneous than those of both DesAb-HSA-D3 and Nb10, and the CDR3 of DesAb-HSA-D3 is slightly more dynamic than that of Nb10, which is the original scaffold used for the design of DesAb-HSA-D3. These differences underline the challenges in the rational design of antibodies by revealing the presence of conformational substates likely to have different binding properties and to generate a high entropic cost upon binding.

Exploring PI3Kγ binding preference with Eganelisib, Duvelisib, and Idelalisib via energetic, pharmacophore and dissociation pathway analyses

Phosphatidylinositol 3-kinase (PI3K) is the central regulator of cellular functions and is suggested as a target for various diseases; thus, effective PI3K inhibitors provide a promising opportunity for the pharmaceutical intervention of many diseases. Among them, PI3Kγ has received more attention because of its essential role in immune signaling. However, the development of novel selective PI3Kγ inhibitors is a major challenge due to the high sequence homology across the class I PI3K isoforms. Therefore, understanding the substrate specificity and receptor-ligand interaction of PI3Kγ would be an appropriate strategy for the rational design of potent γ-selective inhibitors. In this study, by combining various molecular modeling approaches (including classic and enhanced sampling molecular dynamics (MD) simulations, end-point binding free energy calculations, and pharmacophore models), three quinolinone core-containing inhibitors, Idelalisib/CAL-101, Duvelisib/IPI-145, and Eganelisib/IPI-549, were employed to reveal the selective binding mechanisms targeting PI3Kγ. The classic MD and free energy calculations highlight the significant interaction and some key residues for the selective binding against PI3Kγ. Furthermore, the dissociation pathway analysis based on umbrella sampling simulations reveals that hydrophobic interactions are dominant for binding of the three ligands during the dissociation processes, and cooperation between the P-loop and the ligands always exists in the binding/dissociation process. Finally, the pharmacophore model revealed that IPI-549 contains a unique hydrophobic feature, and PI3Kγ exhibits an important hydrogen bond donor feature of hydrogen amide. These findings may provide some important information for the rational design and optimization of PI3Kγ-selective inhibitors.

Enhancing the Inhomogeneous Photodynamics of Canonical Bacteriophytochrome.

The ability of phytochromes to act as photoswitches in plants and microorganisms depends on interactions between a bilin-like chromophore and a host protein. The interconversion occurs between the spectrally distinct red (Pr) and far-red (Pfr) conformers. This conformational change is triggered by the photoisomerization of the chromophore D-ring pyrrole. In this study, as a representative example of a phytochrome-bilin system, we consider biliverdin IXα (BV) bound to bacteriophytochrome (BphP) from Deinococcus radiodurans. In the absence of light, we use an enhanced sampling molecular dynamics (MD) method to overcome the photoisomerization energy barrier. We find that the calculated free energy (FE) barriers between essential metastable states agree with spectroscopic results. We show that the enhanced dynamics of the BV chromophore in BphP contributes to triggering nanometer-scale conformational movements that propagate by two experimentally determined signal transduction pathways. Most importantly, we describe how the metastable states enable a thermal transition known as the dark reversion between Pfr and Pr, through a previously unknown intermediate state of Pfr. We present the heterogeneity of temperature-dependent Pfr states at the atomistic level. This work paves a way toward understanding the complete mechanism of the photoisomerization of a bilin-like chromophore in phytochromes.

PathDetect-SOM: A Neural Network Approach for the Identification of Pathways in Ligand Binding Simulations.

Understanding the process of ligand–protein recognition is important to unveil biological mechanisms and to guide drug discovery and design. Enhanced-sampling molecular dynamics is now routinely used to simulate the ligand binding process, resulting in the need for suitable tools for the analysis of large data sets of binding events. Here, we designed, implemented, and tested PathDetect-SOM, a tool based on self-organizing maps to build concise visual models of the ligand binding pathways sampled along single simulations or replicas. The tool performs a geometric clustering of the trajectories and traces the pathways over an easily interpretable 2D map and, using an approximate transition matrix, it can build a graph model of concurrent pathways. The tool was tested on three study cases representing different types of problems and simulation techniques. A clear reconstruction of the sampled pathways was derived in all cases, and useful information on the energetic features of the processes was recovered. The tool is available at https://github.com/MottaStefano/PathDetect-SOM.