Multicanonical Molecular Dynamics Research Articles

Binding Mechanism between Platelet Glycoprotein and Cyclic Peptide Elucidated by McMD-Based Dynamic Docking.

The cyclic peptide OS1 (amino acid sequence: CTERMALHNLC), which has a disulfide bond between both termini cysteine residues, inhibits complex formation between the platelet glycoprotein Ibα (GPIbα) and the von Willebrand factor (vWF) by forming a complex with GPIbα. To study the binding mechanism between GPIbα and OS1 and, therefore, the inhibition mechanism of the protein-protein GPIbα-vWF complex, we have applied our multicanonical molecular dynamics (McMD)-based dynamic docking protocol starting from the unbound state of the peptide. Our simulations have reproduced the experimental complex structure, although the top-ranking structure was an intermediary one, where the peptide was bound in the same location as in the experimental structure; however, the β-switch of GPIbα attained a different conformation. Our analysis showed that subsequent refolding of the β-switch results in a more stable binding configuration, although the transition to the native configuration appears to take some time, during which OS1 could dissociate. Our results show that conformational changes in the β-switch are crucial for successful binding of OS1. Furthermore, we identified several allosteric binding sites of GPIbα that might also interfere with vWF binding, and optimization of the peptide to target these allosteric sites might lead to a more effective inhibitor, as these are not dependent on the β-switch conformation.

Enhanced Coarse-Grained Molecular Dynamics Simulation with a Smoothed Hybrid Potential Using a Neural Network Model.

In all-atom (AA) molecular dynamics (MD) simulations, the rugged energy profile of the force field makes it challenging to reproduce spontaneous structural changes in biomolecules within a reasonable calculation time. Existing coarse-grained (CG) models, in which the energy profile is set to a global minimum around the initial structure, are unsuitable to explore the structural dynamics between metastable states far away from the initial structure without any bias. In this study, we developed a new hybrid potential composed of an artificial intelligence (AI) potential and minimal CG potential related to the statistical bond length and excluded volume interactions to accelerate the transition dynamics while maintaining the protein character. The AI potential is trained by energy matching using a diverse structural ensemble sampled via multicanonical (Mc) MD simulation and the corresponding AA force field energy, profile of which is smoothed by energy minimization. By applying the new methodology to chignolin and TrpCage, we showed that the AI potential can predict the AA energy with significantly high accuracy, as indicated by a correlation coefficient (R-value) between the true and predicted energies exceeding 0.89. In addition, we successfully demonstrated that CGMD simulation based on the smoothed hybrid potential can significantly enhance the transition dynamics between various metastable states while preserving protein properties compared to those obtained with conventional CGMD and AAMD.

Elucidation of binding mechanism, affinity and complex structure between mWT1 tumor-associated antigen peptide and HLA-A*24:02.

We have applied our advanced computational and experimental methodologies to investigate the complex structure and binding mechanism of a modified Wilms' Tumor 1 (mWT1) protein epitope to the understudied Asian-dominant allele HLA-A*24:02 (HLA-A24) in aqueous solution. We have applied our developed multicanonical molecular dynamics (McMD)-based dynamic docking method to analyze the binding pathway and mechanism, which we verified by comparing the highest probability structures from simulation with our experimentally solved X-ray crystal structure. Subsequent path sampling MD simulations elucidated the atomic details of the binding process and indicated that first an encounter complex is formed between the N-terminal's positive charge of the 9-residue mWT1 fragment peptide and a cluster of negative residues on the surface of HLA-A24, with the MHC molecule preferring a predominantly closed conformation. The peptide first binds to this closed MHC conformation, forming an encounter complex, after which the binding site opens due to increased entropy of the binding site, allowing the peptide to bind to form the native complex structure. Further sequence and structure analyses also suggest that although the peptide loading complex would help with stabilizing the MHC molecule, the binding depends in a large part on the intrinsic affinity between the MHC molecule and the antigen peptide. Finally, our computational tools and analyses can be of great benefit to study the binding mechanism of different MHC types to their antigens, where it could also be useful in the development of higher affinity variant peptides and for personalized medicine. This article is protected by copyright. All rights reserved.

N-Terminal-Driven Binding Mechanism of an Antigen Peptide to Human Leukocyte Antigen-A*2402 Elucidated by Multicanonical Molecular Dynamic-Based Dynamic Docking and Path Sampling Simulations.

We have applied our advanced multicanonical molecular dynamics (McMD)-based dynamic docking methodology to investigate the binding mechanism of an HIV-1 Nef protein epitope to the Asian-dominant allele human leukocyte antigen (HLA)-A*2402. Even though pMHC complex formation [between a Major histocompatibility complex (MHC) class I molecule, which is encoded by an HLA allele, and an antigen peptide] is one of the fundamental processes of the adaptive human immune response, its binding mechanism has not yet been well studied, partially due to the high allelic variation of HLAs in the population. We have used our developed McMD-based dynamic docking method and have successfully reproduced the native complex structure, which is located near the free energy global minimum. Subsequent path sampling MD simulations elucidated the atomic details of the binding process and indicated that the peptide binding is initially driven by the highly positively charged N-terminus of the peptide that is attracted to the various negatively charged residues on the MHC molecule's surface. Upon nearing the pocket, the second tyrosine residue of the peptide anchors the peptide by strongly binding to the B-site of the MHC molecule via hydrophobic driven interactions, resulting in a very strong bound complex structure. Our methodology can be effectively used to predict the bound complex structures between MHC molecules and their antigens to study their binding mechanism in close detail, which would help with the development of new vaccines against cancers, as well as viral infections such as HIV and COVID-19.

Cyclosporin A: Conformational Complexity and Chameleonicity.

The chameleonic behavior of cyclosporin A (CsA) was investigated through conformational ensembles employing multicanonical molecular dynamics simulations that could sample the cis and trans isomers of N-methylated amino acids; these assessments were conducted in explicit water, dimethyl sulfoxide, acetonitrile, methanol, chloroform, cyclohexane (CHX), and n-hexane (HEX) using AMBER ff03, AMBER10:EHT, AMBER12:EHT, and AMBER14:EHT force fields. The conformational details were discussed employing the free-energy landscapes (FELs) at T = 300 K; it was observed that the experimentally determined structures of CsA were only a part of the conformational space. Comparing the ROESY measurements in CHX-d12 and HEX-d14, the major conformations in those apolar solvents were essentially the same as that in CDCl3 except for the observation of some sidechain rotamers. The effects of the metal ions on the conformations, including the cis/trans isomerization, were also investigated. Based on the analysis of FELs, it was concluded that the AMBER ff03 force field best described the experimentally derived conformations, indicating that CsA intrinsically formed membrane-permeable conformations and that the metal ions might be the key to the cis/trans isomerization of N-methylated amino acids before binding a partner protein.

Accurate Binding Configuration Prediction of a G-Protein-Coupled Receptor to Its Antagonist Using Multicanonical Molecular Dynamics-Based Dynamic Docking.

We have performed dynamic docking between a prototypic G-protein-coupled receptor (GPCR) system, the β2-adrenergic receptor, and its antagonist, alprenolol, using one of the enhanced conformation sampling methods, multicanonical molecular dynamics (McMD), which does not rely on any prior knowledge for the definition of the reaction coordinate. Although we have previously applied our McMD-based dynamic docking protocol to various globular protein systems, its application to GPCR systems would be difficult because of their complicated design, which include a lipid bilayer, and because of the difficulty in sampling the configurational space of a binding site that exists deep inside the GPCR. Our simulations sampled a wide array of ligand-bound and ligand-unbound structures, and we measured 427 binding events during our 48 μs production run. Analysis of the ensemble revealed several stable and meta-stable structures, where the most stable structure at the global free energy minimum matches the experimental one. Additional canonical MD simulations were used for refinement and validation of the structures, revealing that most of the intermediates are sufficiently stable to trap the ligand in these intermediary states and furthermore validated our prediction results. Given the difficulty in reaching the orthosteric binding site, chemical optimization of the compound for the second ranking configuration, which binds near the pocket's entrance, might lead to a high-affinity allosteric inhibitor. Accordingly, we show that the application of our methodology can be used to provide crucial insights for the rational design of drugs that target GPCRs.

Cryptic-site binding mechanism of medium-sized Bcl-xL inhibiting compounds elucidated by McMD-based dynamic docking simulations

We have performed multicanonical molecular dynamics (McMD) based dynamic docking simulations to study and compare the binding mechanism between two medium-sized inhibitors (ABT-737 and WEHI-539) that bind to the cryptic site of Bcl-xL, by exhaustively sampling the conformational and configurational space. Cryptic sites are binding pockets that are transiently formed in the apo state or are induced upon ligand binding. Bcl-xL, a pro-survival protein involved in cancer progression, is known to have a cryptic site, whereby the shape of the pocket depends on which ligand is bound to it. Starting from the apo-structure, we have performed two independent McMD-based dynamic docking simulations for each ligand, and were able to obtain near-native complex structures in both cases. In addition, we have also studied their interactions along their respective binding pathways by using path sampling simulations, which showed that the ligands form stable binding configurations via predominantly hydrophobic interactions. Although the protein started from the apo state, both ligands modulated the pocket in different ways, shifting the conformational preference of the sub-pockets of Bcl-xL. We demonstrate that McMD-based dynamic docking is a powerful tool that can be effectively used to study binding mechanisms involving a cryptic site, where ligand binding requires a large conformational change in the protein to occur.

Dynamic Docking Using Multicanonical Molecular Dynamics: Simulating Complex Formation at the Atomistic Level.

Multicanonical molecular dynamics (McMD)-based dynamic docking has been applied to predict the native binding configurations for several protein receptors and their ligands. Due to the enhanced sampling capabilities of McMD, it can exhaustively sample bound and unbound ligand configurations, as well as receptor conformations, and thus enables efficient sampling of the conformational and configurational space, not possible using canonical MD simulations. As McMD samples a wide configurational space, extensive analysis is required to study the diverse ensemble consisting of bound and unbound structures. By projecting the reweighted ensemble onto the first two principal axes obtained via principal component analysis of the multicanonical ensemble, the free energy landscape (FEL) can be obtained. Further analysis produces representative structures positioned at the local minima of the FEL, where these structures are then ranked by their free energy. In this chapter, we describe our dynamic docking methodology, which has successfully reproduced the native binding configuration for small compounds, medium-sized compounds, and peptide molecules.

Exhaustive search of the configurational space of heat-shock protein 90 with its inhibitor by multicanonical molecular dynamics based dynamic docking.

Multicanonical molecular dynamics based dynamic docking was used to exhaustively search the configurational space of an inhibitor binding to the N-terminal domain of heat-shock protein 90 (Hsp90). The obtained structures at 300 K cover a wide structural ensemble, with the top two clusters ranked by their free energy coinciding with the native binding site. The representative structure of the most stable cluster reproduced the experimental binding configuration, but an interesting conformational change in Hsp90 could be observed. The combined effects of solvation and ligand binding shift the equilibrium from a preferred loop-in conformation in the unbound state to an α-helical one in the bound state for the flexible lid region of Hsp90. Thus, our dynamic docking method is effective at predicting the native binding site while exhaustively sampling a wide configurational space, modulating the protein structure upon binding.

Flexible Docking between Enzyme and Its Inhibitor Using Multicanonical MD Simulations and Binding Free Energy Calculations

We have performed dynamic docking simulations between a flexible receptor and a highly flexible ligand by employing multicanonical molecular dynamics (McMD) [1-2]. We have applied our method to predict the native binding configuration and sample the intermediary binding structures between the enzyme β-secretase 1 with a wide binding pocket and its medium-sized inhibitor 3MR [3]. Representative structures located at free energy minima obtained from McMD were taken and subjected to canonical simulations to refine and validate them, reproducing the native complex structure in high agreement with the experimental data. In addition, the binding free energy was estimated by umbrella sampling (US) simulations along representative pathways obtained from the McMD ensemble, followed by weighted histogram analysis to estimate the affinity, which also reproduced the experimental inhibitory affinity. The sampled ensemble by the US simulations smoothly connected the bound and unbound states, refining the binding pathway while staying true to the McMD ensemble. Interestingly, the loss of interactions between the two molecules along the pathway were clearly shown in the free energy landscape, reiterating the fundamental importance of atomistic interactions to the binding affinity between receptor and drug. [1] N. Kamiya et al. Proteins 70, 41-53 (2008). [2] G.-J. Bekker et al. J. Chem. Theory Comput. 13, 2389-2399 (2017). [3] G.-J. Bekker et al. J. Phys. Chem. B 123, 2479-2490 (2019).

MyPresto/omegagene 2020: a molecular dynamics simulation engine for virtual-system coupled sampling.

The molecular dynamics (MD) method is a promising approach for investigating the molecular mechanisms of microscopic phenomena. In particular, generalized ensemble MD methods can efficiently explore the conformational space with a rugged free-energy surface. However, the implementation and acquisition of technical knowledge for each generalized ensemble MD method are not straightforward for end-users. Here, we present a new version of the myPresto/omegagene software, which is an MD simulation engine tailored for a series of generalized ensemble methods, which are virtual-system coupled multicanonical MD (V-McMD), virtual-system coupled adaptive umbrella sampling (V-AUS), and virtual-system coupled canonical MD (VcMD). This program has been applied in several studies analyzing free-energy landscapes of a variety of molecular systems with all-atom simulations. The updated version provides new functionality for coarse-grained simulations powered by the hydrophobicity scale method. The software package includes a step-by-step tutorial document for enhanced conformational sampling of the poly-glutamine (poly-Q) oligomer expressed as a one-bead per residue model. The myPresto/omegagene software is freely available at the following URL: https://github.com/kotakasahara/omegagene under the Apache2 license.

Effects of number of parallel runs and frequency of bias-strength replacement in generalized ensemble molecular dynamics simulations

Background The generalized ensemble approach with the molecular dynamics (MD) method has been widely utilized. This approach usually has two features. (i) A bias potential, whose strength is replaced during a simulation, is applied. (ii) Sampling can be performed by many parallel runs of simulations. Although the frequency of the bias-strength replacement and the number of parallel runs can be adjusted, the effects of these settings on the resultant ensemble remain unclear. Method In this study, we performed multicanonical MD simulations for a foldable mini-protein (Trp-cage) and two unstructured peptides (8- and 20-residue poly-glutamic acids) with various settings. Results As a result, running many short simulations yielded robust results for the Trp-cage model. Regarding the frequency of the bias-potential replacement, although using a high frequency enhanced the traversals in the potential energy space, it did not promote conformational changes in all the systems.

Conformation and Permeability: Cyclic Hexapeptide Diastereomers.

Conformational ensembles of eight cyclic hexapeptide diastereomers in explicit cyclohexane, chloroform, and water were analyzed by multicanonical molecular dynamics (McMD) simulations. Free-energy landscapes (FELs) for each compound and solvent were obtained from the molecular shapes and principal component analysis at T = 300 K; detailed analysis of the conformational ensembles and flexibility of the FELs revealed that permeable compounds have different structural profiles even for a single stereoisomeric change. The average solvent-accessible surface area (SASA) in cyclohexane showed excellent correlation with the cell permeability, whereas this correlation was weaker in chloroform. The average SASA in water correlated with the aqueous solubility. The average polar surface area did not correlate with cell permeability in these solvents. A possible strategy for designing permeable cyclic peptides from FELs obtained from McMD simulations is proposed.

Multimodal Structural Distribution of the p53 C-Terminal Domain upon Binding to S100B via a Generalized Ensemble Method: From Disorder to Extradisorder.

Intrinsically disordered regions (IDRs) of a protein employ a flexible binding manner when recognizing a partner molecule. Moreover, it is recognized that binding of IDRs to a partner molecule is accompanied by folding, with a variety of bound conformations often being allowed in formation of the complex. In this study, we investigated a fragment of the disordered p53 C-terminal domain (CTDf) that interacts with one of its partner molecules, S100B, as a representative IDR. Although the 3D structure of CTDf in complex with S100B has been previously reported, the specific interactions remained controversial. To clarify these interactions, we performed generalized ensemble molecular dynamics (MD) simulations (virtual-system coupled multicanonical MD, termed V-McMD), which enable effective conformational sampling beyond that provided by conventional MD. These simulations generated a multimodal structural distribution for our system including CTDf and S100B, indicating that CTDf forms a variety of complex structures upon binding to S100B. We confirmed that our results are consistent with chemical shift perturbations and nuclear Overhauser effects that were observed in previous studies. Furthermore, we calculated the conformational entropy of CTDf in bound and isolated (free) states. Comparison of these CTDf entropies indicated that the disordered CTDf shows further increase in conformational diversity upon binding to S100B. Such entropy gain by binding may comprise an important feature of complex formation for IDRs.

Dynamic Docking of a Medium-Sized Molecule to Its Receptor by Multicanonical MD Simulations.

A medium-sized and highly flexible inhibitor to the enzyme β-secretase 1 (BACE), which produces the amyloid β-peptide by cleavage of its precursor protein, was dynamically docked into the large and wide catalytic cleft of BACE that binds to the amyloid-precursor by employing multicanonical molecular dynamics (McMD) simulations. We applied our method to predict the native binding configuration and sample the intermediary structures connecting this natively bound state to the unbound one. Representative structures located at free energy minima obtained from McMD were taken and subjected to canonical simulations to refine and validate them, reproducing the native complex structure in agreement with the experimental data in the most stable structure, i.e., the one at the global minimum. In addition, the binding free energy was estimated by umbrella sampling (US) simulations along representative pathways obtained from the McMD ensemble, followed by weighted histogram analysis to estimate the affinity, which also reproduced the experimental inhibitory affinity. Interestingly, the loss of interactions between the two molecules along the pathway was clearly shown in the free energy landscape, reiterating the fundamental importance of atomistic interactions to the binding affinity between receptor and drug compound. The sampled ensemble by the US simulations smoothly connected the bound and unbound states, refining the binding pathway while staying true to the McMD ensemble.

Accurate Prediction of Complex Structure and Affinity for a Flexible Protein Receptor and Its Inhibitor.

In order to predict the accurate binding configuration as well as the binding affinity for a flexible protein receptor and its inhibitor drug, enhanced sampling with multicanonical molecular dynamics (McMD) simulation and thermodynamic integration (TI) were combined as a general drug docking method. CDK2, cyclin-dependent kinase 2, is involved in the cell cycle regulation. Malfunctions in CDK2 can cause tumorigenesis, and thus it is a potential drug target. Here, we performed a long McMD simulation for docking the inhibitor CS3 to CDK2 starting from the unbound structure. Subsequently, a potential binding/unbinding pathway was given from the multicanonical ensemble, and the binding free energy was readily computed by TI along the pathway. Using this combination, the correct binding configuration of CS3 to CDK2 was obtained, and its affinity coincided well with the experimental value.

Variation of free-energy landscape of the p53 C-terminal domain induced by acetylation: Enhanced conformational sampling.

The C‐terminal domain (CTD) of tumor suppressor protein p53 is an intrinsically disordered region that binds to various partner proteins, where lysine of CTD is acetylated/nonacetylated and histidine neutralized/non‐neutralized. Because of the flexibility of the unbound CTD, a free‐energy landscape (FEL) is a useful quantity for determining its statistical properties. We conducted enhanced conformational sampling of CTD in the unbound state via virtual system coupled multicanonical molecular dynamics, in which the lysine was acetylated or nonacetylated and histidine was charged or neutralized. The fragments were expressed by an all‐atom model and were immersed in an explicit solvent. The acetylation and charge‐neutralization varied FEL greatly, which might be convenient to exert a hub property. The acetylation slightly enhanced alpha‐helix structures that are more compact than sheet/loop conformations. The charge‐neutralization produced hairpins. Additionally, circular dichroism experiments confirmed the computational results. We propose possible binding mechanisms of CTD to partners by investigating FEL. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Dynamics of the Extended String-Like Interaction of TFIIE with the p62 Subunit of TFIIH

General transcription factor II E (TFIIE) contains an acid-rich region (residues 378–393) in its α-subunit, comprising 13 acidic and two hydrophobic (Phe387 and Val390) residues. Upon binding to the p62 subunit of TFIIH, the acidic region adopts an extended string-like structure on the basic groove of the pleckstrin homology domain (PHD) of p62, and inserts Phe387 and Val390 into two shallow pockets in the groove. Here, we have examined the dynamics of this interaction by NMR and molecular dynamics (MD) simulations. Although alanine substitution of Phe387 and/or Val390 greatly reduced binding to PHD, the binding mode of the mutants was similar to that of the wild-type, as judged by the chemical-shift changes of the PHD. NMR relaxation dispersion profiles of the interaction exhibited large amplitudes for residues in the C-terminal half-string in the acidic region (Phe387, Glu388, Val390, Ala391, and Asp392), indicating a two-site binding mode: one corresponding to the final complex structure, and one to an off-pathway minor complex. To probe the off-pathway complex structure, an atomically detailed free-energy landscape of the binding mode was computed by all-atom multicanonical MD. The most thermodynamically stable cluster corresponded to the final complex structure. One of the next stable clusters was the off-pathway structure cluster, showing the reversed orientation of the C-terminal half-string on the PHD groove, as compared with the final structure. MD calculations elucidated that the C-terminal half-acidic-string forms encounter complexes mainly around the positive groove region with nearly two different orientations of the string, parallel and antiparallel to the final structure. Interestingly, the most encountered complexes exhibit a parallel-like orientation, suggesting that the string has a tendency to bind around the groove in the proper orientation with the aid of Phe387 and/or Val390 to proceed smoothly to the final complex structure.

Protein-Protein Interactions and Secondary Structure affect Helix Stability in Crowded Environments

The dense, heterogeneous cellular environment is known to affect protein stability through interactions with other biomacromolecules. The effect of excluded volume due to these biomolecules, also known as crowding agents, on a protein of interest, or test protein, has long been known to increase the stability of a test protein. Recently, it has been recognized that attractive protein-crowder interactions play an important role. These interactions affect protein stability and can destabilize the test protein. Here, we use multicanonical molecular dynamics and a coarse-grained protein model to study the folding thermodynamics of a small helical test protein in the presence of crowding agents that are themselves proteins. In order to separate the roles of crowder hydrophobicity and secondary structure, two series of simulations were conducted. Each series covered a range of crowder hydrophobicities: in one series the crowding agents were alpha-helices; in the second, they were beta-hairpins. Our results show that the stability of the test protein depends on both the hydrophobicity and secondary structure of the surrounding biomolecules. For hydrophilic crowding agents, the peptide is stabilized through entropic factors, for moderately hydrophilic crowders, the peptide is destabilized through crowder-induced stabilization of its unfolded states, and for hydrophobic crowders, the peptide is either stabilized through packing of crowding agents around the peptide or destabilized due to beta-sheet formation.

β-Hairpin Crowding Agents Affect α-Helix Stability in Crowded Environments.

The dense, heterogeneous cellular environment is known to affect protein stability. It is now recognized that attractive "quinary" interactions with other biomacromolecules in the cell, referred to as the crowding agents, play a significant role in determining the stability of the protein of interest or test protein. These attractive interactions can reduce or overcome the stabilizing effect of the excluded volume of the crowding agents. However, the roles of specific interactions, such as hydrogen bonding and side chain-side chain hydrophobic interactions, are still unclear. Here, we use molecular simulation to investigate the roles played by hydrophobic interactions and hydrogen bonding between a small helical test protein and equally sized crowding agent proteins in a fixed β-hairpin configuration. The test protein and crowding agents are represented by a coarse-grained protein model, and we use multicanonical molecular dynamics to study the folding thermodynamics of the test protein. Our results confirm that the stability of the test protein depends on the hydrophobicity of the crowding agents and that the stability of the test protein is reduced through favorable side chain-side chain interactions that preferentially stabilize the unfolded states. In addition, we show that when the intermolecular hydrophobic interactions are more favorable than the intramolecular hydrophobic interactions, the β-rich crowding agents can completely destabilize the test protein, causing it to adopt configurations with increased β-content and preventing it from forming its native helical state. Similarities between our results and those seen in the formation of amyloid fibrils are also discussed.