Biased Molecular Dynamics Simulations Research Articles

Searching Peptide Conformational Space.

We have performed a near complete analysis of the conformational space in terms of minima and transition structures for four small peptide models with a force field energy function. There is a clear trend that minima having a large difference in structure, as measured by the distance in torsional space, are rarely connected by a single transition structure. There is a similar trend that activation energies for conformational transitions correlate with structure differences, such that small conformational changes occur with low energy barriers and vice versa. This suggests that a systematic search for low energy conformational transition structures should focus on pairs of minima that are structurally similar. Eigenvectors from diagonalization of force constant matrices at minima are better at describing conformational transitions than vibrational normal modes, as verified both by overlaps with geometry difference vectors and results from biased molecular dynamics simulations.

Use of Umbrella Sampling to Calculate the Entrance/Exit Pathway for Z-Pro-Prolinal Inhibitor in Prolyl Oligopeptidase.

Prolyl oligopeptidase (POP), a member of the prolyl endopeptidase family, is known to play a role in several neurological disorders. Its primary function is to cleave a wide range of small oligopeptides, including neuroactive peptides. We have used force biased molecular dynamics simulation to study the binding mechanism of POP. We examined three possible binding pathways using Steered Molecular Dynamics (SMD) and Umbrella Sampling (US) on a crystal structure of porcine POP with bound Z-pro-prolinal (ZPP). Using SMD, an exit pathway between the first and seventh blade of the β-propeller domain of POP was found to be a nonviable route. US on binding pathways through the β-propeller tunnel and the TYR190-GLN208 flexible loop at the interface between both POP domains allowed us to isolate the flexible loop pathway as the most probable. Further analysis of that pathway suggests a long-range covariation of the interdomain H-bond network, which indicates the possibility of large-scale domain reorientation observed in bacterial homologues and hypothesized to also occur in human POP.

Accelerating chemical reactions: Exploring reactive free-energy surfaces using accelerated ab initio molecular dynamics

A biased potential molecular dynamics simulation approach, accelerated molecular dynamics (AMD), has been implemented in the framework of ab initio molecular dynamics for the study of chemical reactions. Using two examples, the double proton transfer reaction in formic acid dimer and the hypothetical adiabatic ring opening and subsequent rearrangement reactions in methylenecyclopropane, it is demonstrated that ab initio AMD can be readily employed to efficiently explore the reactive potential energy surface, allowing the prediction of chemical reactions and the identification of metastable states. An adaptive variant of the AMD method is developed, which additionally affords an accurate representation of both the free-energy surface and the mechanism associated with the chemical reaction of interest and can also provide an estimate of the reaction rate.

Metadynamics-Based Mechanistic Interpretation of Functional Crosstalk Between Serotonin 2A and Metabotropic Glutamate 2 receptors

Serotonin 2A (2AR) and metabotropic glutamate 2 (mGluR2) receptors have been shown to form a functional and specific heteromeric complex in mammalian brain and in tissue culture preparations with possible implications in the psychotic symptoms of schizophrenia. Unlike non-antipsychotic drugs, clinically effective atypical antipsychotics (e.g., clozapine) that bind to 2AR increase the glutamate-mediated Gi signaling through the 2AR/mGluR2 heterodimer. The molecular mechanisms underlying these allosteric effects and functional crosstalk are unknown. Using molecular dynamics (MD) simulations enhanced with metadynamics, we investigated at the molecular level the conformational changes induced by atypical antipsychotic or non-antipsychotic 2AR ligands in atomistic representations of interacting 2AR and mGluR2 embedded in an explicit lipid-water environment. First, we sampled the conformational transitions from inactive to activated (opsin-like) models of the ligand-free transmembrane regions of 2AR or mGluR2 with adiabatic biased MD simulations. We then reconstructed the free-energy landscape of the 2AR/mGluR2 heterodimer along the pre-determined transition trajectories in the presence of ligands, using a path collective variable approach based on metadynamics. The CHARMM force-field with the CMAP backbone energy correction was used to describe the full systems. All calculations were performed using NAMD enhanced with the Plumed plug-in. Our results suggest that the conformational transitions of 2AR and mGluR2 are populated by different inactive and active metastable states of the receptors, which are differentially stabilized by antipsychotic and non-antipsychotic ligands.

Structural Insights into the Two Sequential Folding Transition States of the PB1 Domain of NBR1 from Φ Value Analysis and Biased Molecular Dynamics Simulations

The PB1 domain of NBR1 folds via a single pathway mechanism involving two sequential energy barriers separated by a high-energy intermediate. The structural ensemble representing each of the two transition states (TS1 and TS2) has been calculated using experimental Φ values and biased molecular dynamics simulations. Both TS1 and TS2 represent compact states (β(TS1) = 0.71, and β(TS2) = 0.93) but are defined by quite different distributions of Φ values, degrees of structural heterogeneity, and nativelike secondary structure. TS1 forms a heterogeneous ensemble of dynamic structures, representing a global collapse of the polypeptide chain around a set of weak nativelike contacts. In contrast, TS2 has a high proportion of nativelike secondary structure, which is reflected in an extensive distribution of high Φ values. Two snapshots along the folding pathway of the PB1 domain reveal insights into the malleability, the solvent accessibility, and the timing of nativelike core packing that stabilizes the folded state.

Downhill binding energy surface of the barnase–barstar complex

We have employed biased molecular dynamics simulations in explicit solvent to characterize the one-dimensional potential of mean force for the dissociation process of the barnase-barstar protein-protein complex. Unbinding of barstar from wild-type barnase was compared with dissociation from four charge-deletion mutants of barnase. Interestingly, we find in all cases that unbinding of barnase and barstar is an uphill process on a smooth, tilted energy landscape. The total free energy difference between the dissociated and bound state was similar for wild-type barnase-barstar and for the R87A mutant of barnase. The values for the three other mutant barnase mutants K27A, R59A, and R83Q were only about half as much. Besides, we have analyzed the conformational dynamics of important residues at the barnase-barstar interface. In the bound state, their conformational fluctuations are reduced relatively to the free state because of the formation of intermolecular contacts. Interestingly, we find that some residues also show decreased mobility at intermediate stages of the unbinding process suggesting that these residues may be involved in the first contacts being formed on binding. © 2010 Wiley Periodicals, Inc. Biopolymers 93: 977-985, 2010.

Putative Active States of a Prototypic G-Protein-Coupled Receptor from Biased Molecular Dynamics

A major current focus of structural work on G-protein-coupled receptors (GPCRs) pertains to the investigation of their active states. However, for virtually all GPCRs, active agonist-bound intermediate states have been difficult to characterize experimentally owing to their higher conformational flexibility, and thus intrinsic instability, as compared to inactive inverse agonist-bound states. In this work, we explored possible activation pathways of the prototypic GPCR bovine rhodopsin by means of biased molecular dynamics simulations. Specifically, we used an explicit atomistic representation of the receptor and its environment, and sampled the conformational transition from the crystal structure of a photoactivated deprotonated state of rhodopsin to the low pH crystal structure of opsin in the presence of 11-trans-retinal, using adiabatic biased molecular dynamics simulations. We then reconstructed the system free-energy landscape along the predetermined transition trajectories using a path collective variable approach based on metadynamics. Our results suggest that the two experimental endpoints of rhodopsin/opsin are connected by at least two different pathways, and that the conformational transition is populated by at least four metastable states of the receptor, characterized by a different amplitude of the outward movement of transmembrane helix 6.

Barnase–Barstar: From first encounter to final complex

Formation of transient protein complexes is an important process in cells. Details of the association process as well as the energy landscapes of association are not well understood. In particular, the nature, height and position of the energy barriers during complexation are debated. Computational studies are well suited for atomistically investigating protein association processes. The Barnase-Barstar complex constitutes a well-studied target for computational studies as a small system with fast association rates. Here, we performed constraint biased Molecular Dynamics simulations along the reaction coordinate reaching from the diffusion regime to the bound complex. We simulated the wild-type and different mutants at different salt concentrations. A structural analysis of our simulation trajectories revealed not a single, but two distinct association patterns dominated by an interplay between two charged contact points near the binding site. Electrostatics and/or mutations influence the relative population of these patterns. Further, we computed the energy landscape of association as PMF (Potential of Mean Force) profiles within a reasonable agreement to experiment. We find a single energy barrier at a distance of ~0.3 nm, which corresponds to the final desolvation transition. Electrostatics has a profound influence on the height of this energy barrier, but not on its position.

A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations.

Trp-cage is a designed 20-residue polypeptide that, in spite of its size, shares several features with larger globular proteins. Although the system has been intensively investigated experimentally and theoretically, its folding mechanism is not yet fully understood. Indeed, some experiments suggest a two-state behavior, while others point to the presence of intermediates. In this work we show that the results of a bias-exchange metadynamics simulation can be used for constructing a detailed thermodynamic and kinetic model of the system. The model, although constructed from a biased simulation, has a quality similar to those extracted from the analysis of long unbiased molecular dynamics trajectories. This is demonstrated by a careful benchmark of the approach on a smaller system, the solvated Ace-Ala3-Nme peptide. For the Trp-cage folding, the model predicts that the relaxation time of 3100 ns observed experimentally is due to the presence of a compact molten globule-like conformation. This state has an occupancy of only 3% at 300 K, but acts as a kinetic trap. Instead, non-compact structures relax to the folded state on the sub-microsecond timescale. The model also predicts the presence of a state at of 4.4 Å from the NMR structure in which the Trp strongly interacts with Pro12. This state can explain the abnormal temperature dependence of the and chemical shifts. The structures of the two most stable misfolded intermediates are in agreement with NMR experiments on the unfolded protein. Our work shows that, using biased molecular dynamics trajectories, it is possible to construct a model describing in detail the Trp-cage folding kinetics and thermodynamics in agreement with experimental data.

Substrate Binding Mechanism of HIV-1 Protease from Explicit-Solvent Atomistic Simulations

The binding mechanism of a peptide substrate (Thr-Ile-Met-Met-Gln-Arg, cleavage site p2-NC of the viral polyprotein) to wild-type HIV-1 protease has been investigated by 1.6 micros biased all-atom molecular dynamics simulations in explicit water. The configuration space has been explored biasing seven reaction coordinates by the bias-exchange metadynamics technique. The structure of the Michaelis complex is obtained starting from the substrate outside the enzyme within a backbone rmsd of 0.9 A. The calculated free energy of binding is -6 kcal/mol, and the kinetic constants for association and dissociation are 1.3 x 10(6) M(-1) s(-1) and 57 s(-1), respectively, consistent with experiments. In the main binding pathway, the flaps of the protease do not open sizably. The substrate slides inside the enzyme cavity from the tight lateral channel. This may contrast with the natural polyprotein substrate which is expected to bind by opening the flaps. Thus, mutations might influence differently the binding kinetics of peptidomimetic ligands and of the natural substrate.

Computational Approach for Ranking Mutant Enzymes According to Catalytic Reaction Rates

A computationally efficient approach for ranking mutant enzymes according to the catalytic reaction rates is presented. This procedure requires the generation and equilibration of the mutant structures, followed by the calculation of partial free energy curves using an empirical valence bond potential in conjunction with biased molecular dynamics simulations and umbrella integration. The individual steps are automated and optimized for computational efficiency. This approach is used to rank a series of 15 dihydrofolate reductase mutants according to the hydride transfer reaction rate. The agreement between the calculated and experimental changes in the free energy barrier upon mutation is encouraging. The computational approach predicts the correct direction of the change in free energy barrier for all mutants, and the correlation coefficient between the calculated and experimental data is 0.82. This general approach for ranking protein designs has implications for protein engineering and drug design.

Exploring the Thermodynamics of Activation Pathways of Bovine Rhodopsin with Fast Molecular Dynamics Simulations

Recent progress made in the understanding of the structure and function of G-protein coupled receptors (GPCRs) points to compelling biophysical and crystallographic evidence for distinct conformational changes of the prototypical family member bovine rhodopsin/opsin upon illumination. To investigate possible activation pathways of bovine rhodopsin from an early inactive intermediate (lumirhodopsin) produced by the cis-trans photoisomerization of its retinal chromophore to an activated ligand-bound opsin-like conformation of the protein, we carried out several independent biased molecular dynamics simulations of the receptor in an explicit dipalmitoylphosphatidylcholine (DPPC) membrane bilayer. The simulations revealed significantly different activation pathways between inactive and opsin-like activated forms of bovine rhodopsin, whose thermodynamic behavior was evaluated by metadynamics using as reaction coordinates both the position along the identified representative pathways, and the distance from them. The results point to three common metastable states that are likely to correspond to a photoactivated deprotonated Meta I-380/Meta IIa, and two different Meta IIb-like intermediates of bovine rhodopsin. Inferences from these molecular models are expected to shed new light into mechanisms of allosteric modulation of rhodopsin function.

COMPUTER SIMULATIONS OF THE INTERACTION OF CIGUATOXIN 3C, BREVENAL AND ent-BREVENAL LADDER POLYETHERS WITH A HOMOLOGY MODEL OF THE VOLTAGE-GATED Kv1.5 POTASSIUM CHANNEL

The interaction of ladder polyethers of marine origin, like ciguatoxin 3C and brevenal, as well as hypothetic ent-brevenal, with the human voltage-gated Kv1.5 potassium ion channel is investigated in this work using homology modeling, automated docking, and energy scoring from molecular dynamics (MD) simulations. A 3D homology model of the pore region of the Kv1.5 channel, previously developed from the 2.9 Å resolution crystal structure of the mammalian Kv1.2 channel — which has a very similar pore sequence — is used here. While ciguatoxin 3C did not enter the pore, both brevenal and ent-brevenal were found into the pore, the latter one with the best score. Binding is attended by notable strain in the ligands, and the corresponding energy increase was evaluated for ent-brevenal by self consistent field (SCF) and density functional theory (DFT) procedures. Egress of ent-brevenal from the pore, as a microscopical reversal of the ingress, was investigated by a smart form of biased MD simulations. While this study indicates ample room and attractive interactions for both brevenal and its enantiomer into the pore, whether these molecules will be found to inhibit voltage-gated potassium ion currents depends upon the barriers in the real system to access the pore, with their thermodynamic and kinetic requirements.

Protein Folding Pathways Revealed by Essential Dynamics Sampling.

The characterization of the protein folding process represents one of the major challenges in molecular biology. Here, a method to simulate the folding process of a protein to its native state is reported, the essential dynamics sampling (EDS) method, and is successfully applied to detecting the correct folding pathways of two small proteins, the all-β SH3 domain of Src tyrosine kinase transforming protein (SH3) and the α/β B1 domain of streptococcal protein G (GB1). The main idea of the method is that a subset of the natural modes of fluctuation in the native state is key in directing the folding process. A biased molecular dynamics simulation is performed, in which the restrained degrees of freedom are chosen among those obtained by a principal component, or essential dynamics, analysis of the positional fluctuations of the Cα atoms in the native state. Successful folding is obtained if the restraints are applied only to the eigenvectors with lowest eigenvalues, representing the most rigid quasi-constraint motions. If the essential eigenvectors, the ones accounting for most of the variance, are used, folding is not successful. These results clearly show that the eigenvectors with lowest eigenvalues contain the main mechanical information necessary to drive the folding process, while the essential eigenvectors represent the large concerted motions which can occur without folding/unfolding the protein.

Atomistic Insights into Rhodopsin Activation from a Dynamic Model

Rhodopsin, the light sensitive receptor responsible for blue-green vision, serves as a prototypical G protein-coupled receptor (GPCR). Upon light absorption, it undergoes a series of conformational changes that lead to the active form, metarhodopsin II (META II), initiating a signaling cascade through binding to the G protein transducin (G(t)). Here, we first develop a structural model of META II by applying experimental distance restraints to the structure of lumi-rhodopsin (LUMI), an earlier intermediate. The restraints are imposed by using a combination of biased molecular dynamics simulations and perturbations to an elastic network model. We characterize the motions of the transmembrane helices in the LUMI-to-META II transition and the rearrangement of interhelical hydrogen bonds. We then simulate rhodopsin activation in a dynamic model to study the path leading from LUMI to our META II model for wild-type rhodopsin and a series of mutants. The simulations show a strong correlation between the transition dynamics and the pharmacological phenotypes of the mutants. These results help identify the molecular mechanisms of activation in both wild type and mutant rhodopsin. While static models can provide insights into the mechanisms of ligand recognition and predict ligand affinity, a dynamic model of activation could be applicable to study the pharmacology of other GPCRs and their ligands, offering a key to predictions of basal activity and ligand efficacy.

Simulation of the Coupling between Nucleotide Binding and Transmembrane Domains in the ATP Binding Cassette Transporter BtuCD

The nucleotide-induced structural rearrangements in ATP binding cassette (ABC) transporters, leading to substrate translocation, are largely unknown. We have modeled nucleotide binding and release in the vitamin B 12 importer BtuCD using perturbed elastic network calculations and biased molecular dynamics simulations. Both models predict that nucleotide release decreases the tilt between the two transmembrane domains and opens the cytoplasmic gate. Nucleotide binding has the opposite effect. The observed coupling may be relevant for all ABC transporters because of the conservation of nucleotide binding domains and the shared role of ATP in ABC transporters. The rearrangements in the cytoplasmic gate region do not provide enough space for B 12 to diffuse from the transporter pore into the cytoplasm, which could suggest that peristaltic forces are needed to exclude B 12 from the transporter pore.

Characterization of the molten globule state of retinol‐binding protein using a molecular dynamics simulation approach

Retinol-binding protein transports retinol, and circulates in the plasma as a macromolecular complex with the protein transthyretin. Under acidic conditions retinol-binding protein undergoes a transition to the molten globule state, and releases the bound retinol ligand. A biased molecular dynamics simulation method has been used to generate models for the ensemble of conformers populated within this molten globule state. Simulation conformers, with a radius of gyration at least 1.1 A greater than that of the native state, contain on average 37%beta-sheet secondary structure. In these conformers the central regions of the two orthogonal beta-sheets that make up the beta-barrel in the native protein are highly persistent. However, there are sizable fluctuations for residues in the outer regions of the beta-sheets, and large variations in side chain packing even in the protein core. Significant conformational changes are seen in the simulation conformers for residues 85-104 (beta-strands E and F and the E-F loop). These changes give an opening of the retinol-binding site. Comparisons with experimental data suggest that the unfolding in this region may provide a mechanism by which the complex of retinol-binding protein and transthyretin dissociates, and retinol is released at the cell surface.

Dissociation of an antiviral compound from the internal pocket of human rhinovirus 14 capsid

WIN antiviral compounds bind human rhinovirus, as well as enterovirus and parechovirus, in an internal cavity located within the viral protein capsid. Access to the buried pocket necessitates deviation from the average viral protein structure identified by crystallography. We investigated the dissociation of WIN 52084 from the pocket in human rhinovirus 14 by using an adiabatic, biased molecular dynamics simulation method. Multiple dissociation trajectories are used to characterize the pathway. WIN 52084 exits between the polypeptide chain near the ends of betaC and betaH in a series of steps. Small, transient packing defects in the protein are sufficient for dissociation. A number of torsion-angle transitions of the antiviral compound are involved, which suggests that flexibility in antiviral compounds is important for binding. It is interesting to note that dissociation is associated with an increase in the conformational fluctuations of residues never in direct contact with WIN 52084 over the course of dissociation. These residues are N-terminal residues in the viral proteins VP3 and VP4 and are located in the interior of the capsid near the icosahedral 5-fold axis. The observed changes in dynamics may be relevant to structural changes associated with virion uncoating and its inhibition by antiviral compounds.

Forces and energetics of hapten-antibody dissociation: a biased molecular dynamics simulation study

The unbinding of fluorescein from the single-chain Fv fragment of the 4D5Flu antibody is investigated by biased molecular dynamics with an implicit solvation model. To obtain statistically meaningful results, a large number of unbinding trajectories are calculated; they involve a total simulation time of more than 200 ns. Simulations are carried out with a time-dependent perturbation and in the presence of a constant force. The two techniques, which provide complementary information, induce unbinding by favoring an increase in the distance between the ligand and the antibody. This distance is an appropriate progress variable for the dissociation reaction and permits direct comparison of the unbinding forces in the simulations with data from atomic force microscopy (AFM). The time-dependent perturbation generates unfolding pathways that are close to equilibrium and can be used to reconstruct the mean force; i.e. the derivative of the potential of mean force, along the reaction coordinate. This is supported by an analysis of the overall unbinding profile and the magnitude of the mean force, which are similar to those of the unbinding force (i.e. the external force due to the time-dependent perturbation) averaged over several unbinding events.The multiple simulations show that unbinding proceeds along a rather wel-defined pathway for a broad range of effective pulling speeds. Initially, there is a distortion of the protein localized in the C-terminal region followed by the fluorescein exit from the binding site. This occurs in steps that involve breaking of specific electrostatic and van der Waals interactions. It appears that the simulations do not explore the same barriers as those measured in the AFM experiments because of the much higher unfolding speed in the former. The dependence of the force on the logarithm of the loading rate is linear and the slope is higher than in the AFM, in agreement with experiment in other systems, where different slopes were observed for different regimes. Based on the unbinding events, mutations in the 4D5Flu antigen binding site are predicted to result in significant changes in the unbinding force.

Analyzing biased Monte Carlo and molecular dynamics simulations

Using the maximum likelihood method, a formalism is derived to analyze a series of biased Monte Carlo or molecular dynamics simulations. The formalism is applied to different examples, in particular the estimation of thermodynamic properties of molecular systems such as potentials of mean force and free energy differences. The formalism is shown to be a generalization of existing methods that are known to be efficient. For the derivation, it is assumed that the values of the distribution function of subsequent states that are analyzed are uncorrelated.