Concept: Molecular dynamics
Post-translational modifications (PTMs) play a key role in numerous cellular processes by directly affecting structure, dynamics and interaction networks of target proteins. Despite their importance, our understanding of protein PTMs at the atomistic level is still largely incomplete. Molecular dynamics (MD) simulations, which provide high-resolution insight into biomolecular function and underlying mechanisms, are in principle ideally suited to tackle this problem. However, because of the challenges associated with the development of novel MD parameters and a general lack of suitable computational tools for incorporating PTMs in target protein structures, MD simulations of post-translationally modified proteins have historically lagged significantly behind the studies of unmodified proteins. Here, we present Vienna-PTM web server (http://vienna-ptm.univie.ac.at), a platform for automated introduction of PTMs of choice to protein 3D structures (PDB files) in a user-friendly visual environment. With 256 different enzymatic and non-enzymatic PTMs available, the server performs geometrically realistic introduction of modifications at sites of interests, as well as subsequent energy minimization. Finally, the server makes available force field parameters and input files needed to run MD simulations of modified proteins within the framework of the widely used GROMOS 54A7 and 45A3 force fields and GROMACS simulation package.
Bcl-X is a member of Bcl-2 family of proteins involved in the regulation of intrinsic pathway of apoptosis. Its overexpression in many human cancers makes it an important target for anti-cancer drugs. Bcl-X interacts with the BH3 domain of several pro-apoptotic Bcl-2 partners. This helical bundle protein has a pronounced hydrophobic groove which acts as a binding region for the BH3 domains. Eight independent molecular dynamics simulations of the apo/holo forms of Bcl-X were carried out to investigate the behavior of solvent-exposed hydrophobic groove. The simulations used either a twin-range cut-off or particle mesh Ewald (PME) scheme to treat long-range interactions. Destabilization of the BH3 domain-containing helix H2 was observed in all four twin-range cut-off simulations. Most of the other major helices remained stable. The unwinding of H2 can be related to the ability of Bcl-X to bind diverse BH3 ligands. The loss of helical character can also be linked to the formation of homo- or hetero-dimers in Bcl-2 proteins. Several experimental studies have suggested that exposure of BH3 domain is a crucial event before they form dimers. Thus unwinding of H2 seems to be functionally very important. The four PME simulations, however, revealed a stable helix H2. It is possible that the H2 unfolding might occur in PME simulations at longer time scales. Hydrophobic residues in the hydrophobic groove are involved in stable interactions among themselves. The solvent accessible surface areas of bulky hydrophobic residues in the groove are significantly buried by the loop LB connecting the helix H2 and subsequent helix. These observations help to understand how the hydrophobic patch in Bcl-X remains stable in the solvent-exposed state. We suggest that both the destabilization of helix H2 and the conformational heterogeneity of loop LB are important factors for binding of diverse ligands in the hydrophobic groove of Bcl-X.
Explicit solvent molecular dynamics simulations have been used to complement preceding experimental and computational studies of folding of guanine quadruplexes (G-DNA). We initiate early stages of unfolding of several G-DNAs by simulating them under no-salt conditions and then try to fold them back using standard excess salt simulations. There is a significant difference between G-DNAs with all-anti parallel stranded stems and those with stems containing mixtures of syn and anti guanosines. The most natural rearrangement for all-anti stems is a vertical mutual slippage of the strands. This leads to stems with reduced numbers of tetrads during unfolding and a reduction of strand slippage during refolding. The presence of syn nucleotides prevents mutual strand slippage; therefore, the antiparallel and hybrid quadruplexes initiate unfolding via separation of the individual strands. The simulations confirm the capability of G-DNA molecules to adopt numerous stable locally and globally misfolded structures. The key point for a proper individual folding attempt appears to be correct prior distribution of syn and anti nucleotides in all four G-strands. The results suggest that at the level of individual molecules, G-DNA folding is an extremely multi-pathway process that is slowed by numerous misfolding arrangements stabilized on highly variable timescales.
Sucrose isomerase NX-5 from Erwiniarhapontici efficiently catalyzes the isomerization of sucrose to isomaltulose (main product) and trehalulose (by-product). To investigate the molecular mechanism controlling sucrose isomer formation, we determined the crystal structures of native NX-5 and its mutant complexes E295Q/sucrose and D241A/glucose at 1.70 Å, 1.70 Å and 2.00 Å, respectively. The overall structure and active site architecture of NX-5 resemble those of other reported sucrose isomerases. Strikingly, the substrate binding mode of NX-5 is also similar to that of trehalulose synthase from Pseudomonasmesoacidophila MX-45 (MutB). Detailed structural analysis revealed the catalytic RXDRX motif and the adjacent 10-residue loop of NX-5 and isomaltulose synthase PalI from Klebsiella sp. LX3 adopt a distinct orientation from those of trehalulose synthases. Mutations of the loop region of NX-5 resulted in significant changes of the product ratio between isomaltulose and trehalulose. The molecular dynamics simulation data supported the product specificity of NX-5 towards isomaltulose and the role of the loop(330-339) in NX-5 catalysis. This work should prove useful for the engineering of sucrose isomerase for industrial carbohydrate biotransformations.
Proper treatment of nonbonded interactions is essential for the accuracy of molecular dynamics (MD) simulations, especially in studies of lipid bilayers. The use of the CHARMM36 force field (C36 FF) in different MD simulation programs can result in disagreements with published simulations performed with CHARMM due to differences in the protocols used to treat the long-range and 1-4 nonbonded interactions. In this study, we systematically test the use of the C36 lipid FF in NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM. A wide range of Lennard-Jones (LJ) cutoff schemes and integrator algorithms were tested to find the optimal simulation protocol to best match bilayer properties of six lipids with varying acyl chain saturation and head groups. MD simulations of a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer were used to obtain the optimal protocol for each program. MD simulations with all programs were found to reasonably match the DPPC bilayer properties (surface area per lipid, chain order parameters, and area compressibility modulus) obtained using the standard protocol used in CHARMM as well as from experiments. The optimal simulation protocol was then applied to the other five lipid simulations and resulted in excellent agreement between results from most simulation programs as well as with experimental data. AMBER compared least favorably with the expected membrane properties, which appears to be due to its use of the hard-truncation in the LJ potential versus a force-based switching function used to smooth the LJ potential as it approaches the cutoff distance. The optimal simulation protocol for each program has been implemented in CHARMM-GUI. This protocol is expected to be applicable to the remainder of the additive C36 FF including the proteins, nucleic acids, carbohydrates, and small molecules.
In this data article we provide topologies and force field parameters files for molecular dynamics simulations of lipids in the OPLS-aa force field using the GROMACS package. This is the first systematic parameterization of lipid molecules in this force field. Topologies are provided for four phosphatidylcholines: saturated DPPC, mono-cis unsaturated POPC and DOPC, and mono-trans unsaturated PEPC. Parameterization of the phosphatidylcholines was achieved in two steps: first, we supplemented the OPLS force field parameters for DPPC with new parameters for torsion angles and van der Waals parameters for the carbon and hydrogen atoms in the acyl chains, as well as new partial atomic charges and parameters for torsion angles in the phosphatidylcholine and glycerol moieties . Next, we derived parameters for the cis and trans double bonds and the neighboring them single bonds . Additionally, we provide GROMACS input files with parameters describing simulation conditions (md.mdp), which are strongly recommended to be used with these lipids models. The data are associated with the research article “Cis and trans unsaturated phosphatidylcholine bilayers: a molecular dynamics simulation study”  and provided as supporting materials.
While most existing theoretical studies on the borophene are based on first-principles calculations, the present work presents molecular dynamics simulations for the lattice dynamical and mechanical properties in borophene. The obtained mechanical quantities are in good agreement with previous first-principles calculations. The key ingredients for these molecular dynamics simulations are the two efficient empirical potentials developed in the present work for the interaction of borophene with low-energy triangular structure. The first one is the valence force field model, which is developed with the assistance of the phonon dispersion of borophene. The valence force field model is a linear potential, so it is rather efficient for the calculation of linear quantities in borophene. The second one is the Stillinger-Weber potential, whose parameters are derived based on the valence force field model. The Stillinger-Weber potential is applicable in molecular dynamics simulations of nonlinear physical or mechanical quantities in borophene.
Molecular dynamics (MD) simulation in the explicit water is performed to study the interaction mechanism of trypsin-ligand binding under the AMBER force field and polarized protein-specific charge (PPC) force field combined the new developed highly efficient interaction entropy (IE) method for calculation of entropy change. And the detailed analysis and comparison of the results of MD simulation for two trypsin-ligand systems show that the root-mean-square deviation (RMSD) of backbone atoms, B-factor, intra-protein and protein-ligand hydrogen bonds are more stable under PPC force field than AMBER force field. Our results demonstrate that the IE method is superior than the traditional normal mode (Nmode) method in the calculation of entropy change and the calculated binding free energy under the PPC force field combined with the IE method is more close to the experimental value than other three combinations (AMBER-Nmode, AMBER-IE and PPC-Nmode). And three critical hydrogen bonds between trypsin and ligand are broken under AMBER force field. However, they are well preserved under PPC force field. Detailed binding interactions of ligands with trypsin are further analyzed. The present work demonstrates that the polarized force field combined the highly efficient IE method is critical in MD simulation and free energy calculation.
Proteins persist longer in the fossil record than DNA, but the longevity, survival mechanisms and substrates remain contested. Here, we demonstrate the role of mineral binding in preserving the protein sequence in ostrich (Struthionidae) eggshell, including from the palaeontological sites of Laetoli (3.8 Ma) and Olduvai Gorge (1.3 Ma) in Tanzania. By tracking protein diagenesis back in time we find consistent patterns of preservation, demonstrating authenticity of the surviving sequences. Molecular dynamics simulations of struthiocalcin-1 and -2, the dominant proteins within the eggshell, reveal that distinct domains bind to the mineral surface. It is the domain with the strongest calculated binding energy to the calcite surface that is selectively preserved. Thermal age calculations demonstrate that the Laetoli and Olduvai peptides are 50 times older than any previously authenticated sequence (equivalent to ~16 Ma at a constant 10°C).
Human immunodeficiency virus type 1 (HIV-1) infection is highly dependent on its capsid. The capsid is a large container, made of ∼1,300 proteins with altogether 4 million atoms. Although the capsid proteins are all identical, they nevertheless arrange themselves into a largely asymmetric structure made of hexamers and pentamers. The large number of degrees of freedom and lack of symmetry pose a challenge to studying the chemical details of the HIV capsid. Simulations of over 64 million atoms for over 1 μs allow us to conduct a comprehensive study of the chemical-physical properties of an empty HIV-1 capsid, including its electrostatics, vibrational and acoustic properties, and the effects of solvent (ions and water) on the capsid. The simulations reveal critical details about the capsid with implications to biological function.