Convergence and sampling efficiency in replica exchange simulations of peptide folding in explicit solvent

Replica exchange methods (REMs) are increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. However, despite having been shown to be very effective on model systems, the application of REM in complex systems such as for the simulation of protein and peptide folding in explicit solvent has not been objectively tested in detail. Here we present a comparison of conventional MD and temperature replica exchange MD (T-REMD) simulations of a beta-heptapeptide in explicit solvent. This system has previously been shown to undergo reversible folding on the time scales accessible to MD simulation and thus allows a direct one-to-one comparison of efficiency. The primary properties compared are the free energy of folding and the relative populations of different conformers as a function of temperature. It is found that to achieve a similar degree of precision T-REMD simulations starting from a random set of initial configurations were approximately an order of magnitude more computationally efficient than a single 800 ns conventional MD simulation for this system at the lowest temperature investigated (275 K). However, whereas it was found that T-REMD simulations are more than four times more efficient than multiple independent MD simulations at one temperature (300 K) the actual increase in conformation sampling was only twofold. The overall gain in efficiency using REMD resulted primarily from the ordering of different conformational states over temperature, as opposed to a large increase of conformational sampling. It is also shown that in this system exchanges are accepted primarily based on (random) fluctuations within the solvent and are not strongly correlated with the instantaneous peptide conformation raising questions in regard to the efficiency of T-REMD in larger systems.     

A molecular dynamics study of the pentacyclo-undecane cage amino acid tripeptide

In this paper, we report on the conformational profile of the pentacyclo-undecane (PCU) cage tripeptide carried out by molecular dynamics (MD) simulation using water as an explicit solvent. The MD solution phase studies carried on the model peptide analogues (A)=Ac–Ala–Ala–Ala–NHMe; (B)=Ac–Cage–Cage–Cage–NHMe; (C)=Ac–Ala–Cage–Ala–NHMe and (D)=Ac–Ala–Pro–Ala–NHMe, are used as a complimentary technique to the corresponding gas phase simulated annealing (SA) study previously carried out in our laboratory. No significant structural changes were observed over the MD trajectories. However, the results reported here provide further evidence that the (PCU) cage amino acid exhibits C_7eq, C_7aq, _R and _L conformations, and the theoretical results suggest that the PCU cage amino acid is a strong β-turn inducer. These results support the prediction that when the PCU cage residues are in the (i) and (i+2) positions, the β-turn can be extended in either direction to form anti-parallel β-pleated sheets, thereby forming the basis of the mechanism for the folding back of the chain in a cross-β-turn structure.     

Formation of left-handed helices in hybrid peptide oligomers with cis beta-sugar amino acid and L-Ala as building blocks

Residue based control of specific helical folding is explored in hybrid peptide oligomers consisting of alternating L-Ala and cis-beta-furanoid sugar amino acid (FSAA) residues as building blocks; two series of these hybrid oligomers are designed, synthesized and extensively characterized by using NMR, CD, FT-IR and MD simulation studies; results show the co-existence of left-handed 11- and 14/15-helical conformations in these short oligomers of Boc-(alpha/beta) and Boc-(beta/alpha) series.     

Folding behaviors of apocytochrome b5 and its mutants: Insights from high temperature molecular dynamics simulations

Apocytochrome b₅ (apocyt b₅) with heme removal from the heme-binding core 1, and its mutants with amino acid replaced in the hydrophobic core 2, namely apocyt b₅ Y7P, P81A and H15R/S20E, have been subjected to molecular dynamics (MD) simulation at high temperature (500 K) for elucidating their folding behaviors. The early events upon thermal induced unfolding were found to be in good agreement with available experimental results, and the lowest stability of Y7P was predicted among the four apoproteins. The influences of these key residues on protein folding behavior were compared directly at an atomic level. At the same time, the influences of non-native interactions of hydrogen bonds and salt-bridges on protein stabilities were analyzed in detail. The insights from current MD simulations are valuable for understanding the apoprotein folding and the holoprotein formation in terms of heme proteins.     

Molecular dynamics simulation of folding of a short helical peptide with many charged residues

A molecular dynamics simulation of the folding of conantokin-T (con-T), a short helical peptide with 5 helical turns of 21 amino acids with 10 charged residues, was carried out to examine folding pathways for this peptide and to predict the folding rate. In the 18 trajectories run at 300 K, 16 trajectories folded, with an averaged folding time of approximately 50 ns. Two trajectories did not fold in up to 200 ns simulation. The folded structure in folded trajectories is in good agreement with experimental structure. An analysis of the trajectories showed that, at the beginning of a few nanoseconds, helix formation started from residues 5-9 with assistance of a hydrophobic clustering involving Tyr5, Met8, and Leu9. The peptide formed a U-shape mainly due to charge-charge interactions between charged residues at the N- and C-terminus segments. In the next approximately 10 ns, several nonnative charge-charge interactions were broken and nonnative Gla10-Lys18 (this denotes a salt bridge between Gal10 and Lys18) and/or Gla10-Lys19 interactions appeared more frequently in this folding step and the peptide became a fishhook J-shape. From this structure, the peptide folded to the folded state in 7 of all 16 folded trajectories in approximately 15 ns. Alternatively, in approximately 30 ns, the con-T went to a conformation in an L-shape with 4 helical turns and a kink at the Arg13 and Gla14 segment in the other 9 trajectories. Con-T in the L-shape then required another approximately 15 ns to fold into the folded state. In addition, in overall folding times, the former 7 trajectories folded faster with the total folding times all shorter than 45 ns, while the latter 9 trajectories folded at a time longer than 45 ns, resulting in an average folding time of approximately 50 ns. Two major folding intermediates found in 2 nonfolded trajectories are stabilized by charge clusters of 5 and 6 charged residues, respectively. With inclusion of friction and solvent-solvent interactions, which were ignored in the present GB/SA solvation model, the folding time obtained above should be multiplied by a factor of 1.25-1.7 according to a previous, similar simulation study. This results in a folding time of 65-105 ns, slightly shorter than the folding time of 127 ns for an alanine-based peptide of the same length. This suggests that the energy barrier of folding for this type of peptides with many charged residues is slightly lower than alanine-based helical peptides by less than 1 kcal/mol.     

Computer simulations study of biomolecules in non-aqueous or cosolvent/water mixture solutions

Pure organic solvents or mixtures with water are very common environments for studying protein and peptide in solution. These milieu conditions are used either for improving the catalytic performance of enzymes or for studying the effect of solvent on the protein stability and hence gaining insight into the protein folding mechanism. The atomic details of these processes are mainly addressed using computer simulation approaches. In particular, Molecular Dynamics simulation represents the most powerful and versatile tool to investigate the details of solvation processes at atomic level. In the last few years, the number of publications peptide and protein simulations in non-natural environments has proliferated. These studies are providing important contributions to shed light on the nature of non-aqueous solvent effects. In this review, the achievements and the future prospects in this field of computational biochemistry are reviewed by summarizing the most important theoretical results published in the last 10 years.     

Comparison of Secondary Structure Formation Using 10 Different Force Fields in Microsecond Molecular Dynamics Simulations

We have compared molecular dynamics (MD) simulations of a β-hairpin forming peptide derived from the protein Nrf2 with 10 biomolecular force fields using trajectories of at least 1 μs. The total simulation time was 37.2 μs. Previous studies have shown that different force fields, water models, simulation methods, and parameters can affect simulation outcomes. The MD simulations were done in explicit solvent with a 16-mer Nrf2 β-hairpin forming peptide using Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, GROMOS96 53a6, CHARMM27, and OPLS-AA/L force fields. The effects of charge-groups, terminal capping, and phosphorylation on the peptide folding were also examined. Despite using identical starting structures and simulation parameters, we observed clear differences among the various force fields and even between replicates using the same force field. Our simulations show that the uncapped peptide folds into a native-like β-hairpin structure at 310 K when Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, or GROMOS96 53a6 were used. The CHARMM27 simulations were able to form native hairpins in some of the elevated temperature simulations, while the OPLS-AA/L simulations did not yield native hairpin structures at any temperatures tested. Simulations that used charge-groups or peptide capping groups were not largely different from their uncapped counterparts with single atom charge-groups. On the other hand, phosphorylation of the threonine residue located at the β-turn significantly affected the hairpin formation. To our knowledge, this is the first study comparing such a large set of force fields with respect to β-hairpin folding. Such a comprehensive comparison will offer useful guidance to others conducting similar types of simulations.     

Implementation of pi-pi interactions in molecular dynamics simulation

No explicit pi-pi interaction term has been incorporated in the conventional molecular dynamics (MD) simulation programs in spite of its significant role in the folding of biomolecules and the clustering of organic chemicals. In this article, we propose a technique to emphasize the effect of pi-pi interactions using a function of energy and implement it into an MD simulation program. Several trial calculations show that the pi-pi incorporated program gives improved results consistent with experimental data on atom geometry and has no unfavorable interference with the conventional computational framework. This indicates an importance of the explicit consideration of pi-pi interactions in MD simulation.     

Examination of the folding of a short alanine-based helical peptide with salt bridges using molecular dynamics simulation

A molecular dynamics simulation of the folding of a short alanine-based helical peptide of 17 residues with three Glu...Lys (i, i + 4) salt bridge pairs, referred to as the AEK17 peptide, was carried out. The simulation gave an estimated simulation folding time of 2.5 ns, shorter than 12 ns for an alanine-based peptide of 16 residues with three Lys residues only, referred to as the AK16 peptide, simulated previously. After folded, the AEK17 peptide had a helical content of 77%, in excellent agreement with the experimentally determined value of 80%. An examination of the folding pathways of AEK17 indicated that the peptide proceeded via three-turn helix conformations more than the helix-turn-helix conformation in the folding pathways. An analysis of interactions indicated that the formation of hydrogen bonds between Lys residue side chains and backbone carbonyls is a major factor in the abundant conformation of the three-turn helix intermediate. The substitution of three Ala with Glu residues reduces the extent of hydrophobic interaction in alanine-based AK peptides with the result that the breaking of the interactions of Lys epsilon-NH3+(side chain)...C=O(backbone) is a major activation action for the AEK17 to achieve a complete fold, in contrast to the AK16 peptide, in which breaking non-native hydrophobic interaction is the rate-determining step.     

Forcefield evaluation and accelerated molecular dynamics simulation of Zn(II) binding to N-terminus of amyloid-β

We report conventional and accelerated molecular dynamics simulation of Zn(II) bound to the N-terminus of amyloid-β. By comparison against NMR data for the experimentally determined binding mode, we find that certain combinations of forcefield and solvent model perform acceptably in describing the size, shape and secondary structure, and that there is no appreciable difference between implicit and explicit solvent models. We therefore used the combination of ff14SB forcefield and GBSA solvent model to compare the result of different binding modes of Zn(II) to the same peptide, using accelerated MD to enhance sampling and comparing the free peptide simulated in the same way. We show that Zn(II) imparts significant rigidity to the peptide, disrupts the secondary structure and pattern of salt bridges seen in the free peptide, and induces closer contact between residues. Free energy surfaces in 1 or 2 dimensions further highlight the effect of metal coordination on peptide’s spatial extent. We also provide evidence that accelerated MD provides improved sampling over conventional MD by visiting as many or more configurations in much shorter simulation times.     

The beta-peptide hairpin in solution: conformational study of a beta-hexapeptide in methanol by NMR spectroscopy and MD simulation.

The structural and thermodynamic properties of a 6-residue beta-peptide that was designed to form a hairpin conformation have been studied by NMR spectroscopy and MD simulation in methanol solution. The predicted hairpin would be characterized by a 10-membered hydrogen-bonded turn involving residues 3 and 4, and two extended antiparallel strands. The interproton distances and backbone torsional dihedral angles derived from the NMR experiments at room temperature are in general terms compatible with the hairpin conformation. Two trajectories of system configurations from 100-ns molecular-dynamics simulations of the peptide in solution at 298 and 340 K have been analyzed. In both simulations reversible folding to the hairpin conformation is observed. Interestingly, there is a significant conformational overlap between the unfolded state of the peptide at each of the temperatures. As already observed in previous studies of peptide folding, the unfolded state is composed of a (relatively) small number of predominant conformers and in this case lacks any type of secondary-structure element. The trajectories provide an excellent ground for the interpretation of the NMR-derived data in terms of ensemble averages and distributions as opposed to single-conformation interpretations. From this perspective, a relative population of the hairpin conformation of 20% to 30% would suffice to explain the NMR-derived data. Surprisingly, however, the ensemble of structures from the simulation at 340 K reproduces more accurately the NMR-derived data than the ensemble from the simulation at 298 K, a question that needs further investigation.     

Structure of sodiated polyglycines

The intrinsic folding of peptides about a sodium ion has been investigated in detail by using infrared multiple photon dissociation (IRMPD) spectroscopy and a combination of theoretical methods. IRMPD spectroscopy was carried out on sodiated polyglycines G(n)-Na(+) (n=2-8), in both the fingerprint and N-H/O-H stretching regions. Interplay between experimental and computational approaches (classical and quantum) enables us to decipher most structural details. The most stable structures of the small peptides up to G(6)-Na(+) maximize metal-peptide interactions with all peptidic C=O groups bound to sodium. In addition, direct interactions between peptide termini are possible for G(6)-Na(+) and larger polyglycines. The increased flexibility of larger peptides leads to more complex folding and internal peptide structuration through γ or β turns. A structural transition is found to occur between G(6)-Na(+) and G(7)-Na(+), leading to a structure with sodium coordination that becomes tri-dimensional for the latter. This transition was confirmed by H/D exchange experiments on G(n)-Na(+) (n=3-8). The most favorable hydrogen-bonding pattern in G(8)-Na(+) involves direct interactions between the peptide termini and opens the way to salt-bridge formation; however, there is only good agreement between experimental and computational data over the entire spectral range for the charge solvation isomer.     

Impact of graphene-based nanomaterials (GBNMs) on the structural and functional conformations of hepcidin peptide

Graphene-based nanomaterials (GBNMs) are widely used in various industrial and biomedical applications. GBNMs of different compositions, size and shapes are being introduced without thorough toxicity evaluation due to the unavailability of regulatory guidelines. Computational toxicity prediction methods are used by regulatory bodies to quickly assess health hazards caused by newer materials. Due to increasing demand of GBNMs in various size and functional groups in industrial and consumer based applications, rapid and reliable computational toxicity assessment methods are urgently needed. In the present work, we investigate the impact of graphene and graphene oxide nanomaterials on the structural conformations of small hepcidin peptide and compare the materials for their structural and conformational changes. Our molecular dynamics simulation studies revealed conformational changes in hepcidin due to its interaction with GBMNs, which results in a loss of its functional properties. Our results indicate that hepcidin peptide undergo severe structural deformations when superimposed on the graphene sheet in comparison to graphene oxide sheet. These observations suggest that graphene is more toxic than a graphene oxide nanosheet of similar area. Overall, this study indicates that computational methods based on structural deformation, using molecular dynamics (MD) simulations, can be used for the early evaluation of toxicity potential of novel nanomaterials.     

Hit-to-Lead Short Peptides against Dengue Type 2 Envelope Protein: Computational and Experimental Investigations

Data from the World Health Organisation show that the global incidence of dengue infection has risen drastically, with an estimated 400 million cases of dengue infection occurring annually. Despite this worrying trend, there is still no therapeutic treatment available. Herein, we investigated short peptide fragments with a varying total number of amino acid residues (peptide fragments) from previously reported dengue virus type 2 (DENV2) peptide-based inhibitors, DN58wt (GDSYIIIGVEPGQLKENWFKKGSSIGQMF), DN58opt (TWWCFYFCRRHHPFWFFYRHN), DS36wt (LITVNPIVTEKDSPVNIEAE), and DS36opt (RHWEQFYFRRRERKFWLFFW), aided by in silico approaches: peptide–protein molecular docking and 100 ns of molecular dynamics (MD) simulation via molecular mechanics using Poisson–Boltzmann surface area (MMPBSA) and molecular mechanics generalised Born surface area (MMGBSA) methods. A library of 11,699 peptide fragments was generated, subjected to in silico calculation, and the candidates with the excellent binding affinity and shown to be stable in the DI-DIII binding pocket of DENV2 envelope (E) protein were determined. Selected peptides were synthesised using conventional Fmoc solid-phase peptide chemistry, purified by RP-HPLC, and characterised using LCMS. In vitro studies followed, to test for the peptides’ toxicity and efficacy in inhibiting the DENV2 growth cycle. Our studies identified the electrostatic interaction (from free energy calculation) to be the driving stabilising force for the E protein–peptide interactions. Five key E protein residues were also identified that had the most interactions with the peptides: (polar) LYS36, ASN37, and ARG350, and (nonpolar) LEU351 and VAL354; these residues might play crucial roles in the effective binding interactions. One of the peptide fragments, DN58opt_8-13 (PFWFFYRH), showed the best inhibitory activity, at about 63% DENV2 plague reduction, compared with no treatment. This correlates well with the in silico studies in which the peptide possessed the lowest binding energy (−9.0 kcal/mol) and was maintained steadily within the binding pocket of DENV2 E protein during the MD simulations. This study demonstrates the use of computational studies to expand research on lead optimisation of antiviral peptides, thus explaining the inhibitory potential of the designed peptides.     

N,N'-linked oligoureas as foldamers: chain length requirements for helix formation in protic solvent investigated by circular dichroism, NMR spectroscopy, and molecular dynamics

N,N'-linked oligoureas with proteinogenic side chains are peptide backbone mimetics belonging to the gamma-peptide lineage. In pyridine, heptamer 4 adopts a stable helical fold reminiscent of the 2.6(14) helical structure proposed for gamma-peptide foldamers. In the present study, we have used a combination of CD and NMR spectroscopies to correlate far-UV chiroptical properties and conformational preferences of oligoureas as a function of chain length from tetramer to nonamer. Both the intensity of the CD spectra and NMR chemical shift differences between alphaCH2 diastereotopic protons experienced a marked increase for oligomers between four and seven residues. No major change in CD spectra occurred between seven and nine residues, thus suggesting that seven residues could be the minimum length required for stabilizing a dominant conformation. Unexpectedly, in-depth NMR conformational investigation of heptamer 4 in CD3OH revealed that the 2.5 helix probably coexists with partially (un)folded conformations and that Z-E urea isomerization occurs, to some degree, along the backbone. Removing unfavorable electrostatic interactions at the amino terminal end of 4 and adding one H-bond acceptor by acylation with alkyl isocyanate (4 --> 7) was found to reinforce the 2.5 helical population. The stability of the 2.5 helical fold in MeOH is further discussed in light of unrestrained molecular dynamics (MD) simulation. Taken together, these new data provide additional insight into the folding propensity of oligoureas in protic solvent and should be of practical value for the design of helical bioactive oligoureas.     

Efficient evaluation of sampling quality of molecular dynamics simulations by clustering of dihedral torsion angles and Sammon mapping

A direct conformational clustering and mapping approach for peptide conformations based on backbone dihedral angles has been developed and applied to compare conformational sampling of Met-enkephalin using two molecular dynamics (MD) methods. Efficient clustering in dihedrals has been achieved by evaluating all combinations resulting from independent clustering of each dihedral angle distribution, thus resolving all conformational substates. In contrast, Cartesian clustering was unable to accurately distinguish between all substates. Projection of clusters on dihedral principal component (PCA) subspaces did not result in efficient separation of highly populated clusters. However, representation in a nonlinear metric by Sammon mapping was able to separate well the 48 highest populated clusters in just two dimensions. In addition, this approach also allowed us to visualize the transition frequencies between clusters efficiently. Significantly, higher transition frequencies between more distinct conformational substates were found for a recently developed biasing-potential replica exchange MD simulation method allowing faster sampling of possible substates compared to conventional MD simulations. Although the number of theoretically possible clusters grows exponentially with peptide length, in practice, the number of clusters is only limited by the sampling size (typically much smaller), and therefore the method is well suited also for large systems. The approach could be useful to rapidly and accurately evaluate conformational sampling during MD simulations, to compare different sampling strategies and eventually to detect kinetic bottlenecks in folding pathways.     

Symplectic molecular dynamics simulations on specially designed parallel computers

We have developed a computer program for molecular dynamics (MD) simulation that implements the Split Integration Symplectic Method (SISM) and is designed to run on specialized parallel computers. The MD integration is performed by the SISM, which analytically treats high-frequency vibrational motion and thus enables the use of longer simulation time steps. The low-frequency motion is treated numerically on specially designed parallel computers, which decreases the computational time of each simulation time step. The combination of these approaches means that less time is required and fewer steps are needed and so enables fast MD simulations. We study the computational performance of MD simulation of molecular systems on specialized computers and provide a comparison to standard personal computers. The combination of the SISM with two specialized parallel computers is an effective way to increase the speed of MD simulations up to 16-fold over a single PC processor.     

Insight into G-DNA structural polymorphism and folding from sequence and loop connectivity through free energy analysis

The lengths of G-tracts and their connecting loop sequences determine G-quadruplex folding and stability. Complete understanding of the sequence-structure relationships remains elusive. Here, single-loop G-quadruplexes were investigated using explicit solvent molecular dynamics (MD) simulations to characterize the effect of loop length, loop sequence, and G-tract length on the folding topologies and stability of G-quadruplexes. Eight loop types, including different variants of lateral, diagonal, and propeller loops, and six different loop sequences [d0 (i.e., no intervening residues in the loop), dT, dT(2), dT(3), dTTA, and dT(4)] were considered through MD simulation and free energy analysis. In most cases the free energetic estimates agree well with the experimental observations. The work also provides new insight into G-quadruplex folding and stability. This includes reporting the observed instability of the left propeller loop, which extends the rules for G-quadruplex folding. We also suggest a plausible explanation why human telomere sequences predominantly form hybrid-I and hybrid-II type structures in K(+) solution. Overall, our calculation results indicate that short loops generally are less stable than longer loops, and we hypothesize that the extreme stability of sequences with very short loops could possibly derive from the formation of parallel multimers. The results suggest that free energy differences, estimated from MD and free energy analysis with current force fields and simulation protocols, are able to complement experiment and to help dissect and explain loop sequence, loop length, and G-tract length and orientation influences on G-quadruplex structure.