Enhanced Sampling Molecular Dynamics Research Articles

The Operando Nature of Isobutene Adsorbed in Zeolite H-SSZ-13 Unraveled by Machine Learning Potentials Beyond DFT Accuracy.

Unraveling the nature of adsorbed olefins in zeolites is crucial to understand numerous zeolite-catalyzed processes. A well-grounded theoretical description critically depends on both an accurate determination of the potential energy surface (PES) and a reliable account of entropic effects at operating conditions. Herein, we propose a transfer learning approach to perform random phase approximation (RPA) quality enhanced sampling molecular dynamics simulations, thereby approaching chemical accuracy on both the determination and exploration of the PES. The proposed methodology is used to investigate isobutene adsorption in H-SSZ-13 as prototypical system to estimate the relative stability of physisorbed olefins, carbenium ions and surface alkoxide species (SAS) in Brönsted-acidic zeolites. We show that the tert-butyl carbenium ion formation is highly endothermic and no entropic stabilization is observed compared to the physisorbed complex within H-SSZ-13. Hence, its predicted concentration and lifetime are negligible, making a direct experimental observation unlikely. Yet, it remains a shallow minimum on the free energy surface over the whole considered temperature range (273--873 K), being therefore a short-lived reaction intermediate rather than a transition state species.

The Operando Nature of Isobutene Adsorbed in Zeolite H‐SSZ‐13 Unraveled by Machine Learning Potentials Beyond DFT Accuracy

Unraveling the nature of adsorbed olefins in zeolites is crucial to understand numerous zeolite‐catalyzed processes. A well‐grounded theoretical description critically depends on both an accurate determination of the potential energy surface (PES) and a reliable account of entropic effects at operating conditions. Herein, we propose a transfer learning approach to perform random phase approximation (RPA) quality enhanced sampling molecular dynamics simulations, thereby approaching chemical accuracy on both the determination and exploration of the PES. The proposed methodology is used to investigate isobutene adsorption in H‐SSZ‐13 as prototypical system to estimate the relative stability of physisorbed olefins, carbenium ions and surface alkoxide species (SAS) in Brönsted‐acidic zeolites. We show that the tert‐butyl carbenium ion formation is highly endothermic and no entropic stabilization is observed compared to the physisorbed complex within H‐SSZ‐13. Hence, its predicted concentration and lifetime are negligible, making a direct experimental observation unlikely. Yet, it remains a shallow minimum on the free energy surface over the whole considered temperature range (273‐‐873 K), being therefore a short‐lived reaction intermediate rather than a transition state species.

CDR L3 Loop Rearrangement Switches Multispecific SPE-7 IgE Antibody From Hapten to Protein Binding.

The monoclonal IgE antibody SPE-7 was originally raised against a 2,4-dinitrophenyl (DNP) target. Through its ability to adopt multiple conformations, the antibody is capable of binding to a diverse range of small haptens and large proteins. The present study examines a dataset of experimentally determined crystal structures of the SPE-7 antibody to gain insight into the mechanisms that contribute to its multispecificity. With the emergence of more and more therapeutic antibodies against a huge repertoire of different targets, our research could be of great interest for future drug development. We are able to discriminate between the different paratope-binding states in the conformational ensembles obtained by enhanced sampling molecular dynamics simulations, and to calculate their transition timescales and state probabilities. Furthermore, we describe the key residues responsible for discriminating between the different binding capacities and identify a tryptophan in a central position of the CDR L3 loop as the residue of greatest interest. The overall dynamics of the paratope appear to be mainly influenced by the CDR L3 and CDR L1 loops.

Mechanism of Dual-Site Recognition in a Classic DNA Aptamer.

Nucleic acid aptamers possess unique advantages in specific recognition. However, the lack of in-depth investigation into their dynamic recognition mechanisms has restricted their rational design and potential applications in fields such as biosensing and targeted therapy. We herein utilized enhanced sampling molecular dynamics to address affinities of adenosine monophosphate (AMP) to the dual binding sites in the DNA aptamer, focusing on the dynamic recognition mechanism and pathways. The present results indicate that in addition to the widely known intermolecular interactions, inequivalence of chemical environments of the two binding sites leads to slightly higher stability of AMP binding to the site proximal to the aptamer terminus. In the presence of two AMPs captured by the two sites, each binding free energy is enhanced. In particular, an additional hydrogen bond of AMP to A10 is introduced in the dual-site binding complex, which increases the binding energy from -4.25 ± 0.47 to -9.48 ± 0.33 kcal mol-1 in the site close to the loop. For the dual-site recognition process, the free energy landscape and minimum free energy pathway calculations elucidate the crucial role of electrostatic interactions between the AMP phosphate groups and Na+ ions in positively cooperative binding mechanisms.

Permeation mechanisms of hydrogen peroxide and water through Plasma Membrane Intrinsic Protein aquaporins.

Hydrogen peroxide (H2O2) transport by aquaporins (AQP) is a critical feature for cellular redox signaling. However, the H2O2 permeation mechanism through these channels remains poorly understood. Through functional assays, two Plasma membrane Intrinsic Protein (PIP) AQP from Medicago truncatula, MtPIP2;2 and MtPIP2;3 have been identified as pH-gated channels capable of facilitating the permeation of both water (H2O) and H2O2. Employing a combination of unbiased and enhanced sampling molecular dynamics simulations, we investigated the key barriers and translocation mechanisms governing H2O2 permeation through these AQP in both open and closed conformational states. Our findings reveal that both H2O and H2O2 encounter their primary permeation barrier within the selectivity filter (SF) region of MtPIP2;3. In addition to the SF barrier, a second energetic barrier at the NPA (asparagine-proline-alanine) region that is more restrictive for the passage of H2O2 than for H2O, was found. This behavior can be attributed to a dissimilar geometric arrangement and hydrogen bonding profile between both molecules in this area. Collectively, these findings suggest mechanistic heterogeneity in H2O and H2O2 permeation through PIPs.

Dimerization Promotes PKR Activation by Modulating Energetics of αC Helix Conversion between Active and Inactive Conformations.

Protein kinase R (PKR) functions in the eukaryotic innate immune system as a first-line defense against viral infections. PKR binds viral dsRNA, leading to autophosphorylation and activation. In its active state, PKR can phosphorylate its primary substrate, eIF2α, which blocks the initiation of translation in the infected cell. It has been established that PKR activation occurs when the kinase domain dimerizes in a back-to-back configuration. However, the mechanism by which dimerization leads to enzymatic activation is not fully understood. Here, we investigate the structural mechanistic basis and energy landscape for PKR activation, with a focus on the αC helix─a kinase activation and signal integration hub─using all-atom equilibrium and enhanced sampling molecular dynamics simulations. By employing window-exchange umbrella sampling, we compute free-energy profiles of activation, which show that back-to-back dimerization stabilizes a catalytically competent conformation of PKR. Key hydrophobic residues in the homodimer interface contribute to stabilization of the αC helix in an active conformation and the position of its critical glutamate residue. Using linear mutual information analysis, we analyze allosteric communication connecting the protomers' N-lobes and the αC helix dimer interface with the αC helix.

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.

Structural visualization of inhibitor binding in prolyl oligopeptidase.

The association and dissociation of proteins and ligands are crucial in biophysics for potential drug development [Baron and McCammon, Annu. Rev. Phys. Chem. 64, 151-175 (2013)]. However, identifying and characterizing the reaction pathways for these rare events has been a long-standing challenge. Molecular dynamics (MD) simulations are limited in exploring biophysical processes on experimental timescales, so ligand transport processes through complex transient tunnels formed by proteins during dynamics are often simulated using enhanced sampling MD [Rydzewski and Nowak, Phys. Life Rev. 22-23, 58-74 (2017)]. Erroneously identified ligand binding pathways can affect thermodynamic and kinetic characteristics calculated from MD trajectories. A system that has the potential to be a therapeutic target for neurodegenerative diseases is prolyl oligopeptidase (PREP). This is due to its involvement in promoting protein aggregation and disrupting cellular function through affecting protein-protein interactions (PPI). The recent discovery of a new type of PREP inhibitor that targets PPI raises important questions about the diversity of ligand binding pathways in PREP and their impact on protein dynamics [Pätsi et al., J. Med. Chem. 67, 5421-5436 (2024); Kilpeläinen et al., J. Med. Chem. 66, 7475-7496 (2023); and Walczewska-Szewc et al., Phys. Chem. Chem. Phys. 24, 4366-4373 (2022)]. In this article, using results from enhanced sampling MD, we visually present how the binding process in PREP depends on subtle changes in inhibitors, which could be crucial in treating neurodegenerative disorders.

Folding of N-terminally acetylated α-synuclein upon interaction with lipid membranes

α-Synuclein (α-syn) is an abundant presynaptic neuronal protein whose aggregation is strongly associated with Parkinson’s disease. It has been proposed that lipid membranes significantly affect α-syn’s aggregation process. Extensive studies have been conducted to understand the interactions between α-syn and lipid membranes and have demonstrated that the N-terminus plays a critical role. However, the dynamics of the interactions and the conformational transitions of the N-terminus of α-syn at the atomistic scale details are still highly desired. In this study, we performed extensive enhanced sampling molecular dynamics simulations to quantify the folding and interactions of wild-type and N-terminally acetylated α-syn when interacting with lipid structures. We found that N-terminal acetylation significantly increases the helicity of the first few residues in solution or when interacting with lipid membranes. The observations in simulations showed that the binding of α-syn with lipid membranes mainly follows the induced-fit model, where the disordered α-syn binds with the lipid membrane through the electrostatic interactions and hydrophobic contacts with the packing defects; after stable insertion, N-terminal acetylation promotes the helical folding of the N-terminus to enhance the anchoring, thus increasing the binding affinity. We have shown the critical role of the first N-terminal residue methionine for recognition and anchoring to the negatively charged membrane. Although N-terminal acetylation neutralizes the positive charge of Met1 that may affect the electrostatic interactions of α-syn with membranes, the increase in helicity of the N-terminus should compensate for the binding affinity. This study provides detailed insight into the folding dynamics of α-syn’s N-terminus with or without acetylation in solution and upon interaction with lipids, which clarifies how the N-terminal acetylation regulates the affinity of α-syn binding to lipid membranes. It also shows how packing defects and electrostatic effects coregulate the N-terminus of α-syn folding and its interaction with membranes.

Exploring an Intracellular Allosteric Site of CC-Chemokine Receptor4 from 3D Models, Probe Simulations, and Mutagenesis.

We applied our previously developed probe confined dynamic mapping protocol, which combines enhanced sampling molecular dynamics (MD) simulations and fragment-based approaches, to identify the binding site of GSK2239633A (N-[[3-[[3-[(5-chlorothiophen-2-yl)sulfonylamino]-4-methoxyindazol-1-yl]methyl]phenyl]methyl]-2-hydroxy-2-methylpropanamide), a selective CC-chemokine receptor type 4 (CCR4) negative allosteric modulator, using CCR4 homology and AlphaFold models. By comparing the performance across five computational models, we identified conserved (K3108.49 and Y3047.53) and non-conserved (M2436.36) residue hotspots for GSK2239633A binding, which were validated by mutagenesis and bioluminescence resonance energy transfer assay. Further analysis of 3D models and MD simulations highlighted the pair of residues 6.36 and 7.56 that might account for antagonist selectivity among chemokine receptors. Our in silico protocol provides a promising approach for characterizing ligand binding sites in membrane proteins, considering receptor dynamics and adaptability and guiding protein template selection for ligand design.

Dilated cardiomyopathy mutation in beta-cardiac myosin enhances actin activation of the power stroke and phosphate release.

Inherited mutations in human beta-cardiac myosin (M2β) can lead to severe forms of heart failure. The E525K mutation in M2β is associated with dilated cardiomyopathy (DCM) and was found to stabilize the interacting heads motif (IHM) and autoinhibited super-relaxed (SRX) state in dimeric heavy meromyosin. However, in monomeric M2β subfragment 1 (S1) we found that E525K enhances (threefold) the maximum steady-state actin-activated ATPase activity (k cat) and decreases (eightfold) the actin concentration at which ATPase is one-half maximal (K ATPase). We also found a twofold to fourfold increase in the actin-activated power stroke and phosphate release rate constants at 30 μM actin, which overall enhanced the duty ratio threefold. Loaded motility assays revealed that the enhanced intrinsic motor activity translates to increased ensemble force in M2β S1. Glutamate 525, located near the actin binding region in the so-called activation loop, is highly conserved and predicted to form a salt bridge with another conserved residue (lysine 484) in the relay helix. Enhanced sampling molecular dynamics simulations predict that the charge reversal mutation disrupts the E525-K484 salt bridge, inducing conformations with a more flexible relay helix and a wide phosphate release tunnel. Our results highlight a highly conserved allosteric pathway associated with actin activation of the power stroke and phosphate release and suggest an important feature of the autoinhibited IHM is to prevent this region of myosin from interacting with actin. The ability of the E525K mutation to stabilize the IHM likely overrides the enhanced intrinsic motor properties, which may be key to triggering DCM pathogenesis.

Automated collective variable discovery for MFSD2A transporter from molecular dynamics simulations

Biomolecules often exhibit complex free energy landscapes in which long-lived metastable states are separated by large energy barriers. Overcoming these barriers to robustly sample transitions between the metastable states with classical molecular dynamics (MD) simulations presents a challenge. To circumvent this issue, collective variable (CV)-based enhanced sampling MD approaches are often employed. Traditional CV selection relies on intuition and prior knowledge of the system. This approach introduces bias, which can lead to incomplete mechanistic insights. Thus, automated CV detection is desired to gain a deeper understanding of the system/process. Analysis of MD data with various machine-learning algorithms, such as principal component analysis (PCA), support vector machine, and linear discriminant analysis (LDA) based approaches have been implemented for automated CV detection. However, their performance has not been systematically evaluated on structurally and mechanistically complex biological systems. Here, we applied these methods to MD simulations of the MFSD2A (Major Facilitator Superfamily Domain 2A) lysolipid transporter in multiple functionally relevant metastable states with the goal of identifying optimal CVs that would structurally discriminate these states. Specific emphasis was on the automated detection and interpretive power of LDA-based CVs. We found that LDA methods, which included a novel gradient descent-based multiclass harmonic variant, termed GDHLDA, we developed here, outperform PCA in class separation, exhibiting remarkable consistency in extracting CVs critical for distinguishing metastable states. Furthermore, the identified CVs included features previously associated with conformational transitions in MFSD2A. Specifically, conformational shifts in transmembrane helix 7 and in residue Y294 on this helix emerged as critical features discriminating the metastable states in MFSD2A. This highlights the effectiveness of LDA-based approaches in automatically extracting from MD trajectories CVs of functional relevance that can be used to drive biased MD simulations to efficiently sample conformational transitions in the molecular system.

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.

Emle-engine: A Flexible Electrostatic Machine Learning Embedding Package for Multiscale Molecular Dynamics Simulations.

We present in this work the emle-engine package (https://github.com/chemle/emle-engine)─the implementation of a new machine learning embedding scheme for hybrid machine learning potential/molecular-mechanics (ML/MM) dynamics simulations. The package is based on an embedding scheme that uses a physics-based model of the electronic density and induction with a handful of tunable parameters derived from in vacuo properties of the subsystem to be embedded. This scheme is completely independent of the in vacuo potential and requires only the positions of the atoms of the machine learning subsystem and the positions and partial charges of the molecular mechanics environment. These characteristics allow emle-engine to be employed in existing QM/MM software. We demonstrate that the implemented electrostatic machine learning embedding scheme (named EMLE) is stable in enhanced sampling molecular dynamics simulations. Through the calculation of free energy surfaces of alanine dipeptide in water with two different ML options for the in vacuo potential and three embedding models, we test the performance of EMLE. When compared to the reference DFT/MM surface, the EMLE embedding is clearly superior to the MM one based on fixed partial charges. The configurational dependence of the electronic density and the inclusion of the induction energy introduced by the EMLE model leads to a systematic reduction in the average error of the free energy surface when compared to MM embedding. By enabling the usage of EMLE embedding in practical ML/MM simulations, emle-engine will make it possible to accurately model systems and processes that feature significant variations in the charge distribution of the ML subsystem and/or the interacting environment.

Enhanced Sampling Molecular Dynamics Simulations Reveal Transport Mechanism of Glycoconjugate Drugs through GLUT1.

Glucose transporters GLUT1 belong to the major facilitator superfamily and are essential to human glucose uptake. The overexpression of GLUT1 in tumor cells designates it as a pivotal target for glycoconjugate anticancer drugs. However, the interaction mechanism of glycoconjugate drugs with GLUT1 remains largely unknown. Here, we employed all-atom molecular dynamics simulations, coupled to steered and umbrella sampling techniques, to examine the thermodynamics governing the transport of glucose and two glycoconjugate drugs (i.e., 6-D-glucose-conjugated methane sulfonate and 6-D-glucose chlorambucil) by GLUT1. We characterized the specific interactions between GLUT1 and substrates at different transport stages, including substrate recognition, transport, and releasing, and identified the key residues involved in these procedures. Importantly, our results described, for the first time, the free energy profiles of GLUT1-transporting glycoconjugate drugs, and demonstrated that H160 and W388 served as important gates to regulate their transport via GLUT1. These findings provide novel atomic-scale insights for understanding the transport mechanism of GLUT1, facilitating the discovery and rational design of GLUT1-targeted anticancer drugs.

Gas adsorption and framework flexibility of CALF-20 explored via experiments and simulations

In 2021, Svante, in collaboration with BASF, reported successful scale up of CALF-20 production, a stable MOF with high capacity for post-combustion CO2 capture which exhibits remarkable stability towards water. CALF-20’s success story in the MOF commercialisation space provides new thinking about appropriate structural and adsorptive metrics important for CO2 capture. Here, we combine atomistic-level simulations with experiments to study adsorptive properties of CALF-20 and shed light on its flexible crystal structure. We compare measured and predicted CO2 and water adsorption isotherms and explain the role of water-framework interactions and hydrogen bonding networks in CALF-20’s hydrophobic behaviour. Furthermore, regular and enhanced sampling molecular dynamics simulations are performed with both density-functional theory (DFT) and machine learning potentials (MLPs) trained to DFT energies and forces. From these simulations, the effects of adsorption-induced flexibility in CALF-20 are uncovered. We envisage this work would encourage development of other MOF materials useful for CO2 capture applications in humid conditions.

Deep learning path-like collective variable for enhanced sampling molecular dynamics.

Several enhanced sampling techniques rely on the definition of collective variables to effectively explore free energy landscapes. The existing variables that describe the progression along a reactive pathway offer an elegant solution but face a number of limitations. In this paper, we address these challenges by introducing a new path-like collective variable called the "deep-locally non-linear-embedding," which is inspired by principles of the locally linear embedding technique and is trained on a reactive trajectory. The variable mimics the ideal reaction coordinate by automatically generating a non-linear combination of features through a differentiable generalized autoencoder that combines a neural network with a continuous k-nearest neighbor selection. Among the key advantages of this method is its capability to automatically choose the metric for searching neighbors and to learn the path from state A to state B without the need to handpick landmarks a priori. We demonstrate the effectiveness of DeepLNE by showing that the progression along the path variable closely approximates the ideal reaction coordinate in toy models, such as the Müller-Brown potential and alanine dipeptide. Then, we use it in the molecular dynamics simulations of an RNA tetraloop, where we highlight its capability to accelerate transitions and estimate the free energy of folding.

Intramolecular and Water Mediated Tautomerism of Solvated Glycine.

Understanding tautomerism and characterizing solvent effects on the dynamic processes pose significant challenges. Using enhanced-sampling molecular dynamics based on state-of-the-art deep learning potentials, we investigated the tautomeric equilibria of glycine in water. We observed that the tautomerism between neutral and zwitterionic glycine can occur through both intramolecular and intermolecular proton transfers. The latter proceeds involving a contact anionic-glycine-hydronium ion pair or separate cationic-glycine-hydroxide ion pair. These pathways with comparable barriers contribute almost equally to the reaction flux.

Base pair dynamics, electrostatics, and thermodynamics at the LTR-III quadruplex:duplex junction

G-quadruplexes (GQs) play key regulatory roles within the human genome and have also been identified to play similar roles in other eukaryotes, bacteria, archaea, and viruses. Human immunodeficiency virus 1, the etiological agent of acquired immunodeficiency syndrome, can form two GQs in its long terminal repeat (LTR) promoter region, each of which act to regulate viral gene expression in opposing manners. The major LTR GQ, called LTR-III, is a distinct hybrid GQ containing a 12-nucleotide duplex loop attached to the quadruplex motif. The resulting quadruplex:duplex junction (QDJ) has been hypothesized to serve as a selective drug targeting site. To better understand the dynamics of this QDJ, we performed conventional and enhanced-sampling molecular dynamics simulations using the Drude-2017 force field. We observed unbiased and reversible formation of additional base pairs in the QDJ, between Ade4:Thy14 and Gua3:Thy14. Both base pairs were electrostatically favored, but geometric constraints within the junction may drive the formation of, and preference for, the Ade4:Thy14 base pair. Finally, we demonstrated that the base pairs are separated only by small energy barriers that may enable transitions between both base-paired states. Together, these simulations provide new insights into the dynamics, electrostatics, and thermodynamics of the LTR-III QDJ.

Comprehensive Analysis of Coupled Proline Cis-Trans States in Bradykinin Using ωBP-REMD Simulations.

It is well-known that proline (Pro) cis-trans isomerization plays a decisive role in the folding and stabilization of proteins. The conformational coupling between isomerization states of different Pro residues in proteins during conformational adaptation processes is not well understood. In the present work, we investigate the coupled cis-trans isomerization of three Pro residues using bradykinin (BK), a partially unstructured nonapeptide hormone, as a model system. We use a recently developed enhanced-sampling molecular dynamics method (ω-bias potential replica exchange molecular dynamics; ωBP-REMD) that allows us to exhaustively sample all combinations of Pro isomer states and obtain converged probability densities of all eight state combinations within 885 ns ωBP-REMD simulations. In agreement with experiment, the all-trans state is seen to be the preferred isomer of zwitterionic aqueous BK. In about a third of its structures, this state presents the characteristic C-terminal β-turn conformation; however, other isomer combinations also contribute significantly to the structural ensemble. Unbiased probabilities can be projected onto the peptide bond dihedral angles of the three Pro residues. This unveils the interdependence of the individual Pro isomerization states, i.e., a possible coupling of the different Pro isomers. The cis/trans equilibrium of a Pro residue can change by up to 2.5 kcal·mol-1, depending on the isomerization state of other Pro residues. For example, for Pro7, the simulations indicate that its cis state becomes favored compared to its trans state when Pro2 is switched from the trans state to the cis state. Our findings demonstrate the efficiency of the ωBP-REMD methodology and suggest that the coupling of Pro isomerization states may play an even more decisive role in larger folded proteins subject to more conformational restraints.