Refining the description of peptide backbone conformations improves protein simulations using the GROMOS 53A6 force field

Conformations and Solvation of Synthetic Polymers in Water by Generalized Born Implicit-Solvent Molecular Dynamics Simulations: Stereoisomers of Poly(acrylic acid) and Poly(methacrylic acid)

We present the GB-OBC model as an approach for implicit-solvent MD simulations of a synthetic macromolecule in water. The model is tested and found to be successful in reproducing the chain dimensions and predicting the free energy of solvation of carboxylic acid vinyl polymers. The influence of stereochemistry and the hydrophobic nature of the polymer was investigated as a function of chain length (20 < N < 600) for poly(acrylic acid) (PAA) and poly(methacrylic acid) (PMA). The dimensionless parameters of the GB-OBC model were parameterized to be applicable to PAA and PMA. Scaling relations for chain dimensions obtained using implicit-solvent MD simulations in this study are in good agreement with those from experiments, theory of solvated chains in good solvents, and all-atom MD simulations in explicit water. Results show that ⟨R_g²⟩/NL² is greater for the atactic chain as compared to the isotactic chain, for PAA as well as PMA. ⟨R_g²⟩/NL² values of chains attain constancy in water for N = 200, with the values being greater for PMA. The PMA chain is conformationally more perturbed than the PAA chain, for both isotactic and atactic stereochemistry. The solvation free energy ΔG_hyd of PAA and PMA in water is negative for all chain lengths (N = 20-600) and becomes more favorable with an increase in molecular weight. The ΔG_hyd values for isotactic and atactic chains are identical at lower values of N but differ slightly for N > 300. Irrespective of the hydrophobic nature of the polymer, the atactic chain is thermodynamically more soluble in water as compared to the isotactic chain. The isotactic chain is less hydrophilic as compared to the atactic chain due to the closer proximity of the COOH groups along the backbone. This implicit solvent method is an effective way to accurately simulate the configurational properties and solvation of synthetic polymers in water.

Comparison of the United- and All-Atom Representations of (Halo)alkanes Based on Two Condensed-Phase Force Fields Optimized against the Same Experimental Data Set

The level of accuracy that can be achieved by a force field is influenced by choices made in the interaction-function representation and in the relevant simulation parameters. These choices, referred to here as functional-form variants (FFVs), include for example the model resolution, the charge-derivation procedure, the van der Waals combination rules, the cutoff distance, and the treatment of the long-range interactions. Ideally, assessing the effect of a given FFV on the intrinsic accuracy of the force-field representation requires that only the specific FFV is changed and that this change is performed at an optimal level of parametrization, a requirement that may prove extremely challenging to achieve in practice. Here, we present a first attempt at such a comparison for one specific FFV, namely the choice of a united-atom (UA) versus an all-atom (AA) resolution in a force field for saturated acyclic (halo)alkanes. Two force-field versions (UA vs AA) are optimized in an automated way using the CombiFF approach against 961 experimental values for the pure-liquid densities ρ_liq and vaporization enthalpies ΔH_vap of 591 compounds. For the AA force field, the torsional and third-neighbor Lennard-Jones parameters are also refined based on quantum-mechanical rotational-energy profiles. The comparison between the UA and AA resolutions is also extended to properties that have not been included as parameterization targets, namely the surface-tension coefficient γ, the isothermal compressibility κ_T, the isobaric thermal-expansion coefficient α_P, the isobaric heat capacity c_P, the static relative dielectric permittivity ϵ, the self-diffusion coefficient D, the shear viscosity η, the hydration free energy ΔG_wat, and the free energy of solvation ΔG_che in cyclohexane. For the target properties ρ_liq and ΔH_vap, the UA and AA resolutions reach very similar levels of accuracy after optimization. For the nine other properties, the AA representation leads to more accurate results in terms of η; comparably accurate results in terms of γ, κ_T, α_P, ϵ, D, and ΔG_che; and less accurate results in terms of c_P and ΔG_wat. This work also represents a first step toward the calibration of a GROMOS-compatible force field at the AA resolution.

Molecular Insight into the Effect of a Single-Nucleotide Polymorphic Variation on the Structure and Dynamics of Methionine Synthase Reductase and Its Association with Neural Tube Defects

Neural tube defects (NTDs) are among the common and severe congenital malformations in neonates. According to a WHO report, nearly three lakh babies are affected per year worldwide by NTDs. Most studies revealed that folate deficiency is the key element to promote NTD with other oligogenic and multifactorial elements. This folate is metabolized by the FOCM (folate one-carbon metabolism) pathway. The most important step in the FOCM pathway is the conversion of methionine to homocysteine, which is guided by the enzyme MTRR. Several single-nucleotide polymorphisms (SNPs) in the MTRR gene are strongly associated with the progression of NTD. A nonsynonymous allelic variant (rs1532268) of the protein leads to a missense mutation at the 202nd position from serine to leucine (S202L) and is associated with a higher disease prevalence in different populations. In our study, this SNP indicates a 2-fold increase in the risk of disease progression (p-value of 0.03; OR 2.76; 95% CI 1.08-7.11). Here, extensive molecular dynamics simulations and interaction network analysis reveal that the change of 202nd serine to leucine alters the structures of the FAD and NAD binding domains, which restricts the ligand binding. The S202L variation alters the functional dynamics that might impede the electron transport chain along the NADP(H)→ FAD→ FMN pathway and hamper phosphorylation by kinases like GSK-3 and CaM-II during the posttranscriptional modification of the protein. The present study provides functional insights into the effect of the genetic variations of the MTRR gene on the NTD disease pathogenesis.

Molecular dynamics simulations of an engineered T4 lysozyme exclude helix to sheet transition, and provide insights into long distance, intra-protein switchable motion

An engineered variant of T4 lysozyme serves as a model for studying induced remote conformational changes in a full protein context. The design involves a duplicated surface helix, flanked by two loops, that switches between two different conformations spanning about 20 Å. Molecular dynamics simulations of the engineered protein, up to 1 μs, rule out α-helix to β-sheet transitions within the duplicated helix as suggested by others. These simulations highlight how the use of different force fields can lead to radical differences in the structure of the protein. In addition, Markov state modeling and transition path theory were employed to map a 6.6 μs simulation for possible early intermediate states and to provide insights into the onset of the switching motion. The putative intermediates involve the folding of one helical turn in the C-terminal loop through energy driven, sequential rearrangement of nearby salt bridges around the key residue Arg63. These results provide a first step towards understanding the energetics and dynamics of a rather complicated intra-protein motion.

Effects of Tyr555 and Trp678 on the processivity of cellobiohydrolase A from Ruminiclostridium thermocellum: A simulation study

Cellobiohydrolase A from Ruminiclostridium thermocellum (Cbh9A) is a processive exoglucanase from family 9 and is an important cellobiohydrolase that hydrolyzes cello-oligosaccharide into cellobiose. Residues Tyr555 and Trp678 considerably affect catalytic activity, but their mechanisms are still unknown. To investigate how the Tyr555 and Trp678 affect the processivity of Cbh9A, conventional molecular dynamics, steered molecular dynamics, and free energy calculation were performed to simulate the processive process of wild type (WT)-Cbh9A, Y555S mutant, and W678G mutant. Analysis of simulation results suggests that the binding free energies between the substrate and WT-Cbh9A are lower than those of Y555S and W678G mutants. The pull forces and energy barrier in Y555S and W678G mutants also reduced significantly during the steered molecular dynamics (SMD) simulation compared with that of the WT-Cbh9A. And the potential mean force calculations showed that the pulling energy barrier of Y555S and W678G mutants is much lower than that of WT-Cbh9A.

Trypsin Binding with Copper Ions Scavenges Superoxide: Molecular Dynamics-Based Mechanism Investigation

Trypsin is a serine protease, which has been proved to be a novel superoxide scavenger. The burst of superoxide induced by polychlorinated biphenyls can be impeded by trypsin in both wild type and sod knockout mutants of Escherichia coli. The experimental results demonstrated that the activities of superoxide scavenging of trypsin were significantly accelerated by Cu ions. Also, with the addition of Cu ions, a new β-sheet (β7) transited from a random coil in the Cu(II)-trypsin (TP) system, which was favorable for the formation of more contacts with other sheets of trypsin. Residue–residue network analysis and the porcupine plots proved that the Cu ion in trypsin strengthened some native interactions among residues, which ultimately resulted in much greater stability of the Cu(II)-TP system. Moreover, compact and stable trypsin structures with Cu ions might be responsible for significantly provoking the activity of superoxide scavenging.

A coupled two-dimensional main chain torsional potential for protein dynamics: generation and implementation

Main chain torsions of alanine dipeptide are parameterized into coupled 2-dimensional Fourier expansions based on quantum mechanical (QM) calculations at M06 2X/aug-cc-pvtz//HF/6-31G** level. Solvation effect is considered by employing polarizable continuum model. Utilization of the M06 2X functional leads to precise potential energy surface that is comparable to or even better than MP2 level, but with much less computational demand. Parameterization of the 2D expansions is against the full main chain torsion space instead of just a few low energy conformations. This procedure is similar to that for the development of AMBER03 force field, except unique weighting factor was assigned to all the grid points. To avoid inconsistency between quantum mechanical calculations and molecular modeling, the model peptide is further optimized at molecular mechanics level with main chain dihedral angles fixed before the calculation of the conformational energy on molecular mechanical level at each grid point, during which generalized Born model is employed. Difference in solvation models at quantum mechanics and molecular mechanics levels makes this parameterization procedure less straightforward. All force field parameters other than main chain torsions are taken from existing AMBER force field. With this new main chain torsion terms, we have studied the main chain dihedral distributions of ALA dipeptide and pentapeptide in aqueous solution. The results demonstrate that 2D main chain torsion is effective in delineating the energy variation associated with rotations along main chain dihedrals. This work is an implication for the necessity of more accurate description of main chain torsions in the future development of ab initio force field and it also raises a challenge to the development of quantum mechanical methods, especially the quantum mechanical solvation models.

Directional Dependence of Hydrogen Bonds: a Density-based Energy Decomposition Analysis and Its Implications on Force Field Development

One well-known shortcoming of widely-used biomolecular force fields is the description of the directional dependence of hydrogen bonding (HB). Here we aim to better understand the origin of this difficulty and thus provide some guidance for further force field development. Our theoretical approaches center on a novel density-based energy decomposition analysis (DEDA) method [J. Chem. Phys., 131, 164112 (2009)], in which the frozen density energy is variationally determined through constrained search. This unique and most significant feature of DEDA enables us to find that the frozen density interaction term is the key factor in determining the HB orientation, while the sum of polarization and charge-transfer components shows very little HB directional dependence. This new insight suggests that the difficulty for current non-polarizable force fields to describe the HB directional dependence is not due to the lack of explicit polarization or charge-transfer terms. Using the DEDA results as reference, we further demonstrate that the main failure coming from the atomic point charge model can be overcome largely by introducing extra charge sites or higher order multipole moments. Among all the electrostatic models explored, the smeared charge distributed multipole model (up to quadrupole), which also takes account of charge penetration effects, gives the best agreement with the corresponding DEDA results. Meanwhile, our results indicate that the van der Waals interaction term needs to be further improved to better model directional hydrogen bonding.

Molecular insight into conformational transition of amyloid β-peptide 42 inhibited by (-)-epigallocatechin-3-gallate probed by molecular simulations

Considerable experimental evidence indicates that (-)-epigallocatechin-3-gallate (EGCG) inhibits the fibrillogenesis of Aβ(42) and alleviates its associated cytotoxicity. However, the molecular mechanism of the inhibition effect of EGCG on the conformational transition of Aβ(42) remains unclear due to the limitations of current experimental techniques. In this work, molecular dynamics simulations and molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) analysis were coupled to better understand the issue. It was found that the direct interactions between EGCG and the peptide are the origin of its inhibition effects. Specifically, EGCG molecules expel water from the surface of the Aβ(42), cluster with each other, and interact directly with the peptide. The results of free energy decomposition calculated by MM-PBSA indicate that the nonpolar term contributes more than 71% to the binding free energy of the EGCG-Aβ(42) complex, while polar interactions (i.e., hydrogen bonding) play a minor role. It was identified that there are 12 important residues of Aβ(42) that strongly interact with EGCG (Phe4, Arg5, Phe19, Phe20, Glu22, Lys28, Gly29, Leu34-Gly37, and Ile41), while nonpolar interactions are mainly provided by the side chains of some hydrophobic residues (Phe, Met and Ile) and the main chains of some nonhydrophobic residues (Lys28 and Gly29). On the contrary, polar interactions are mainly formed by the main chain of Aβ(42), of which the main chains of Gly29 and Gly37 contribute greatly. The work has thus elucidated the molecular mechanism of the inhibition effect of EGCG on the conformational transition of Aβ(42), and the findings are considered critical for exploring more effective agents for the inhibition of Aβ(42) fibrillogenesis.

Using molecular dynamics to study liquid phase behavior: simulations of the ternary sodium laurate/sodium oleate/water system

The prediction of surfactant phase behavior has applications in a wide range of areas. An accurate modeling of liquid phase behavior can aid our understanding of colloidal process or be used to design phases that respond in a defined way to their environment. In this work, we use molecular dynamics to model the phase behavior of the ternary sodium laurate/sodium oleate/water system and compare the simulation results to experimental data. Simulations were performed with the GROMOS 53A6 united-atom force field and cover the entire ternary phase diagram, producing micellar, hexagonal, and lamellar phases. The aggregate simulation time for the 33 simulations performed during this study is 4.4 μs. We find that the simulations were able to model the experimentally observed liquid phase behavior accurately, showing that the carboxylate and lipid parameters of the 53A6 force field give very good quality results for the in silico prediction of liquid system phase behavior.

Conformational state-specific free energy differences by one-step perturbation: protein secondary structure preferences of the GROMOS 43A1 and 53A6 force fields

The one-step free energy perturbation approach can be applied to obtain conformational state-specific free energy differences (FEDs) associated with changes in force field parameters, and thus offers the possibility to consider conformational equilibria during force field parameterization. In this work, using the alanine decapeptide in explicit water solution as a model, the α-helical and β-hairpin state-specific FEDs associated with force field changes between two widely used parameter sets of the GROMOS force field, namely, 43A1 and 53A6, were determined using one-step perturbation. The results mostly deviated by only 1 kJ mol(−1) in absolute or a few percent in relative values from thermodynamic integration results, suggesting that the convergence ranges of one-step perturbation were large enough to cover the substantial changes in nonbonded parameters between the two parameter sets. It was also found that one-step perturbation may give larger errors when the changes from the reference state include a large decrease in van der Waals radius, as indicated by the result for the β-hairpin state-specific free energy change going from 53A6 to 43A1. According to the free energy results, the α-helical state of the alanine decapeptide is destabilized by 15 kJ mol(−1), i.e., 1.5 kJ mol(−1) per residue, relative to the β-hairpin state when going from 43A1 to 53A6, in agreement with previous direct simulations in which native α-helices were often found to be unstable in simulations using 53A6, despite that the 53A6 parameters better reproduce a range of thermodynamic properties of small molecular systems. By applying one-step perturbation to analyze the effects of perturbing individual parameters, the differential stabilization of the two secondary structure states can be traced to the changes in van der Waals parameters, especially a van der Waals parameter involved in third-neighbor interactions. This study provides an example of the efficiency of one-step perturbation in force field development, reducing the computational cost by orders of magnitude.

Using one-step perturbation to predict the effect of changing force-field parameters on the simulated folding equilibrium of a beta-peptide in solution

Computer simulation using molecular dynamics is increasingly used to simulate the folding equilibria of peptides and small proteins. Yet, the quality of the obtained results depends largely on the quality of the force field used. This comprises the solute as well as the solvent model and their energetic and entropic compatibility. It is, however, computational very expensive to perform test simulations for each combination of force-field parameters. Here, we use the one-step perturbation technique to predict the change of the free enthalpy of folding of a beta-peptide in methanol solution due to changing a variety of force-field parameters. The results show that changing the solute backbone partial charges affects the folding equilibrium, whereas this is relatively insensitive to changes in the force constants of the torsional energy terms of the force field. Extending the cut-off distance for nonbonded interactions beyond 1.4 nm does not affect the folding equilibrium. The same result is found for a change of the reaction-field permittivity for methanol from 17.7 to 30. The results are not sensitive to the criterion, e.g., atom-positional RMSD or number of hydrogen bonds, that is used to distinguish folded and unfolded conformations. Control simulations with perturbed Hamiltonians followed by backward one-step perturbation indicated that quite large perturbations still yield reliable results. Yet, perturbing all solvent molecules showed where the limitations of the one-step perturbation technique are met. The evaluated methodology constitutes an efficient tool in force-field development for molecular simulation by reducing the number of required separate simulations by orders of magnitude.

Functional roles of a structural element involving Na+-pi interactions in the catalytic site of T1 lipase revealed by molecular dynamics simulations

Interactions between metal ions and pi systems (metal-pi interactions) are known to confer significant stabilization energy. However, in biological systems, few structures with metal-pi coordination have been determined; thus, its roles must still be elucidated. The cation-pi interactions are not correctly described by current molecular mechanics even when using a polarizable force field, and thus they require quantum mechanical calculations for accurate estimation. However, the huge computational costs of the latter methodologies prohibit long-time molecular dynamics (MD) simulations. Accordingly, we developed a novel scheme to obtain an effective potential for calculating the interaction energy with an accuracy comparable to that of advanced ab initio calculations at the CCSD(T) levels, and with computational costs comparable to those of conventional MM calculations. Then, to elucidate the functional roles of the Na(+)-phenylalanine (Phe) complex in the catalytic site of T1 lipase, we performed MD simulations in the presence/absence of the accurate Na(+)-pi interaction energy. A comparison of these MD simulations revealed that a significantly large enthalpy gain in Na(+)-Phe16 substantially stabilizes the catalytic site, whereas a water molecule could not be substituted for Na(+) for sufficient stabilization energy. Thus, the cation-pi interaction in the lipase establishes a remarkably stable core structure by combining a hydrophobic aromatic ring and hydrophilic residues, of which the latter form the catalytic triad, thereby contributing to large structural changes from the complex with ligands to the free form of the lipase. This is the first report to elucidate the detailed functional mechanisms of Na(+)-pi interactions.

Secondary structure propensities in peptide folding simulations: a systematic comparison of molecular mechanics interaction schemes

We present a systematic study directed toward the secondary structure propensity and sampling behavior in peptide folding simulations with eight different molecular dynamics force-field variants in explicit solvent. We report on the combinational result of force field, water model, and electrostatic interaction schemes and compare to available experimental characterization of five studied model peptides in terms of reproduced structure and dynamics. The total simulation time exceeded 18 mus and included simulations that started from both folded and extended conformations. Despite remaining sampling issues, a number of distinct trends in the folding behavior of the peptides emerged. Pronounced differences in the propensity of finding prominent secondary structure motifs in the different applied force fields suggest that problems point in particular to the balance of the relative stabilities of helical and extended conformations.

Structural stability of electrosprayed proteins: temperature and hydration effects

Electrospray ionization is a gentle method for sample delivery, routinely used in gas-phase studies of proteins. It is crucial for structural investigations that the protein structure is preserved, and a good understanding of how structure is affected by the transition to the gas phase is needed for the tuning of experiments to meet that requirement. Small amounts of residual solvent have been shown to protect the protein, but temperature is important too, although it is not well understood how the latter affects structural details. Using molecular dynamics we have simulated four sparingly hydrated globular proteins (Trp-cage; Ctf, a C-terminal fragment of a bacterial ribosomal protein; ubiquitin; and lysozyme) in vacuum starting at temperatures ranging from 225 K to 425 K. For three of the proteins, our simulations show that a water layer corresponding to 3 A preserves the protein structure in vacuum, up to starting temperatures of 425 K. Only Ctf shows minor secondary structural changes at lower starting temperatures. The structural conservation stems mainly from interactions with the surrounding water. Temperature scales in simulations are not directly translatable into experiments, but the wide temperature range in which we find the proteins to be stable is reassuring for the success of future single particle imaging experiments. The water molecules aggregate in clusters and form patterns on the protein surface, maintaining a reproducible hydrogen bonding network. The simulations were performed mainly using OPLS-AA/L, with cross checks using AMBER03 and GROMOS96 53a6. Only minor differences between the results from the three different force fields were observed.