Force field impact and spin-probe modeling in molecular dynamics simulations of spin-labeled T4 lysozyme

Several attempts have been made to compute electron paramagnetic resonance (EPR) spectra of biomolecules, using motional models or simulated trajectories to describe dynamics. Ideally, the simulated trajectories should capture fast (picosecond) snapshots of spin-probe rotations accurately, while being lengthy enough to ensure a proper Fourier integration of the time-domain signal. It is the interplay of the two criteria that poses computational challenges to the method. In this context, an analysis of the spin-probe and protein conformational sampling and equilibration, with different force fields and with explicit solvent, may be a useful attempt. The present work reports a comparative study of the effect of the molecular dynamics (MD) force field on conformational sampling and equilibration in two spin-labeled T4 lysozyme (T4L) variants, N40C and K48C. Ensembles of 10× 3 ns-trajectories per variant and per force field (OPLS/AMBER and AMBER99) are analyzed for a reliable assessment of convergence and sampling. It is found that subtle site-dependent differences in spin-probe rotations and torsions are more readily captured in the AMBER99 trajectories than in the OPLS/AMBER simulations. On the other hand, sampling and equilibration are found to be better with the OPLS/AMBER force field at equal trajectory lengths.Figure: Left panel: The spin-probe R1 ring and the spin-probe Euler angles , and . Middle panel: Illustration of the diffusion in a cone model for the spin-probe motion: snapshots of helix B and of the R1 ring in N40C, taken at 0.3 ps intervals from AMBER trajectory 1. Right panel: The N40C mutant with the spin label (solid mode), solvated in a cubic box.     

Membrane proteins: molecular dynamics simulations

Molecular dynamics simulations of membrane proteins are making rapid progress, because of new high-resolution structures, advances in computer hardware and atomistic simulation algorithms, and the recent introduction of coarse-grained models for membranes and proteins. In addition to several large ion channel simulations, recent studies have explored how individual amino acids interact with the bilayer or snorkel/anchor to the headgroup region, and it has been possible to calculate water/membrane partition free energies. This has resulted in a view of bilayers as being adaptive rather than purely hydrophobic solvents, with important implications, for example, for interaction between lipids and arginines in the charged S4 helix of voltage-gated ion channels. However, several studies indicate that the typical current simulations fall short of exhaustive sampling, and that even simple protein-membrane interactions require at least ca. 1mus to fully sample their dynamics. One new way this is being addressed is coarse-grained models that enable mesoscopic simulations on multi-mus scale. These have been used to model interactions, self-assembly and membrane perturbations induced by proteins. While they cannot replace all-atom simulations, they are a potentially useful technique for initial insertion, placement, and low-resolution refinement.     

Analysis of ligand binding to proteins using molecular dynamics simulations

This work aims to explore theoretically the molecular mechanisms of ligand binding to proteins through the use of molecular dynamics simulations. The binding of sodium dodecyl sulfate (SDS) to cobra cardio toxin A3 (CTX A3) and thiourea (TOU) to lysozyme have been chosen as the two model systems. Data acquisitions were made by Gromacs software. To begin with, the collisions of ligand molecules with every residue of CTX A3 and lysozyme were evaluated. With this information in hand, the average numbers of collisions with each residue was defined and then assessed. Next, a measure of the affinity of a residue, P_i, referred to as conformational factor, toward a ligand molecule was established. Based on the results provided, all site-making residues for CTX A3 and lysozyme were identified. The results are in good agreement with the experimental data. Finally, based on this method, all site-making residues of bovine carbonic anhydrase (BCA) toward the SDS ligand were predicted.     

Application of molecular dynamics simulations to design a dual-purpose oligopeptide linker sequence for fusion proteins

Proteins are often monitored by combining a fluorescent polypeptide tag with the target protein. However, due to the high molecular weight and immunogenicity of such tags, they are not suitable choices for combining with fusion proteins such as immunotoxins. In this study, we designed a polypeptide sequence with a dual role (it acts as both a linker and a fluorescent probe) to use with fusion proteins. Two common fluorescent tag sequences based on tetracysteine were compared to a commonly used rigid linker as well as our proposed dual-purpose sequence. Computational investigations showed that the dual-purpose sequence was structurally stable and may be a good choice to use as both a linker and a fluorescence marker between two moieties in a fusion protein.     

OmpA: Gating and dynamics via molecular dynamics simulations

Outer membrane proteins (OMPs) of Gram-negative bacteria have a variety of functions including passive transport, active transport, catalysis, pathogenesis and signal transduction. Whilst the structures of ∼ 25 OMPs are currently known, there is relatively little known about their dynamics in different environments. The outer membrane protein, OmpA from Escherichia coli has been studied extensively in different environments both experimentally and computationally, and thus provides an ideal test case for the study of the dynamics and environmental interactions of outer membrane proteins. We review molecular dynamics simulations of OmpA and its homologues in a variety of different environments and discuss possible mechanisms of pore gating. The transmembrane domain of E. coli OmpA shows subtle differences in dynamics and interactions between a detergent micelle and a lipid bilayer environment. Simulations of the crystallographic unit cell reveal a micelle-like network of detergent molecules interacting with the protein monomers. Simulation and modelling studies emphasise the role of an electrostatic-switch mechanism in the pore-gating mechanism. Simulation studies have been extended to comparative models of OmpA homologues from Pseudomonas aeruginosa (OprF) and Pasteurella multocida (PmOmpA), the latter model including the periplasmic C-terminal domain.     

Electrostatic screening in molecular dynamics simulations.

The screened Coulombic potential has been shown to describe satisfactorily equilibrium properties like pK shifts, the effects of charged groups on redox potentials and binding constants of metal ions. To test how well the screening of the electrostatic potential describes the dynamical trajectory of a macromolecular system, a series of comparative simulations have been carried out on a protein system which explicitly included water molecules and a system in vacuo. For the system without solvent the results of using (i) the standard potential form were compared with results of (ii) the potential where the Coulomb term was modified by the inclusion of a distance dependent dielectric, epsilon (r), to model the screening effect of bulk water, and (iii) standard potential modified by reducing the charge on ionized residue side chains. All molecular dynamics simulations have been carried out on bovine pancreatic trypsin inhibitor. Comparisons between the resulting trajectories, averaged structures, hydrogen bonding patterns and properties such as solvent accessible surface area and radius of gyration are described. The results show that the dynamical behaviour of the protein calculated with a screened electrostatic term compares more favourably with the time-dependent structural changes of the full system with explicitly included water than the standard vacuum simulation.     

Conformational analysis and molecular dynamics simulations of maltose

Constrained conformational energy minimizations have been used to calculate an adiabatic (Φ, ψ) potential energy surface for the disaccharide β-maltose. The inclusion of molecular flexibility in the conformational energy analysis of the disaccharide was found to significantly lower the barriers to conformational transitions, as has been observed previously for other systems. Several low energy wells were identified on the adiabatic surface which differ in energy by small amounts and with low absolute barriers separating them, indicating the possibility of a non-negligible equilibrium population distribution in each well. If such a distribution of conformations existed in the physical system, the conformation observed by NMR NOE measurements would thus be a “virtual” conformation. Molecular dynamics simulations of the motions of this molecule in vacuum were also conducted and indicate that the rate of relaxation of the molecule to the adiabatic surface may be slower than the typical timescale of conformational fluctuations. This effect is apparently due to an unphysical persistence of hydrogen bond patterns in vacuum which does not occur in aqueous solution. Trajectories undergoing transitions between wells were calculated and the effects of such conformational transitions upon the ensemble mean structure, such as might be observed in an NMR experiment, were demonstrated.     

Predicting thymine dimerization yields from molecular dynamics simulations

It was recently shown that thymine dimers in the all-thymine oligonucleotide (dT)₁₈ are fully formed in <1 ps after ultraviolet excitation. The speed and low quantum yield of this reaction suggest that only a small fraction of the conformers of this structurally disordered oligonucleotide are in a position to react at the instant of photon absorption. In this work, we explore the hypothesis that conventional molecular dynamics simulations can be used to predict the yield of cyclobutane pyrimidine dimers in DNA. Conformations obtained from simulations of thymidylyl-(3′-5′)-thymidine in various cosolvents were classified as dimerizable or nondimerizable depending on the distance between the C5-C6 double bonds of the adjacent thymine bases and the torsion angle between them. The quantum yield of cyclobutane pyrimidine dimer formation was calculated as the number of dimerizable conformations divided by the total number of conformations. The experimental quantum yields measured in the different solvents were satisfactorily reproduced using physically reasonable values for the two parameters. The mean dimerizable structure computed by averaging all of the dimerizable cis-syn conformations is structurally similar to the actual cis-syn dimer. Compared to the canonical B-form TT step, the most important structural property of a dimerizable conformation is its reduced helical twist angle of 22°.     

Protein flexibility using constraints from molecular dynamics simulations

Proteins are held together in the native state by hydrophobic interactions, hydrogen bonds and interactions with the surrounding water, whose strength as well as spatial and temporal distribution affects protein flexibility and hence function. We study these effects using 10 ns molecular dynamics simulations of pure water and of two proteins, the glutamate receptor ligand binding domain and barnase. We find that most of the noncovalent interactions flicker on and off over typically nanoseconds, and so we can obtain good statistics from the molecular dynamics simulations. Based on this information, a topological network of rigid bonds corresponding to a protein structure with covalent and noncovalent bonds is constructed, with account being taken of the influence of the flickering hydrogen bonds. We define the duty cycle for the noncovalent interactions as the percentage of time a given interaction is present, which we use as an input to investigate flexibility/rigidity patterns, in the algorithm FIRST which constructs and analyses topological networks.     

Time-averaged order parameter restraints in molecular dynamics simulations

A method is described that allows experimental \(S^2\) order parameters to be enforced as a time-averaged quantity in molecular dynamics simulations. The two parameters that characterize time-averaged restraining, the memory relaxation time and the weight of the restraining potential energy term in the potential energy function used in the simulation, are systematically investigated based on two model systems, a vector with one end restrained in space and a pentapeptide. For the latter it is shown that the backbone N–H order parameter of individual residues can be enforced such that the spatial fluctuations of quantities depending on atomic coordinates are not significantly perturbed. The applicability to realistic systems is illustrated for the B3 domain of protein G in aqueous solution.     

Sequence-dependent DNA deformability studied using molecular dynamics simulations

Proteins recognize specific DNA sequences not only through direct contact between amino acids and bases, but also indirectly based on the sequence-dependent conformation and deformability of the DNA (indirect readout). We used molecular dynamics simulations to analyze the sequence-dependent DNA conformations of all 136 possible tetrameric sequences sandwiched between CGCG sequences. The deformability of dimeric steps obtained by the simulations is consistent with that by the crystal structures. The simulation results further showed that the conformation and deformability of the tetramers can highly depend on the flanking base pairs. The conformations of xATx tetramers show the most rigidity and are not affected by the flanking base pairs and the xYRx show by contrast the greatest flexibility and change their conformations depending on the base pairs at both ends, suggesting tetramers with the same central dimer can show different deformabilities. These results suggest that analysis of dimeric steps alone may overlook some conformational features of DNA and provide insight into the mechanism of indirect readout during protein–DNA recognition. Moreover, the sequence dependence of DNA conformation and deformability may be used to estimate the contribution of indirect readout to the specificity of protein–DNA recognition as well as nucleosome positioning and large-scale behavior of nucleic acids.     

Atomistic modeling of apatite-collagen composites from molecular dynamics simulations extended to hyperspace

The preparation of atomistic models of apatite-collagen composite mimicking enamel at length scales in the range of 1–10 nanometers is outlined. This bio-composite is characterized by a peculiar interplay of the collagen triplehelix and the apatite crystal structure. Structural coherence is however only obtained after drastic rearrangements, namely the depletion of protein-protein hydrogen bonds and the incorporation of calcium triangles which are stabilized by salt-bridges with the collagen molecule. Starting from an isolated collagen triple helix and a single-crystalline apatite structure, a composite model is obtained by gradually merging the two components via an additional (hyperspace) coordinate. This approach allows smooth structural relaxation of both components whilst avoiding singularities in potential energy due to atomic overlap.     

Higher order multipole moments for molecular dynamics simulations

In conventional force fields, the electrostatic potential is represented by atom-centred point charges. This choice is in principle arbitrary, but technically convenient. Point charges can be understood as the first term of multipole expansions, which converge with an increasing number of terms towards the accurate representation of the molecular potential given by the electron density distribution. The use of multipole expansions can therefore improve the force field accuracy. Technically, the implementation of atomic multipoles is more involved than the use of point charges. Important points to consider are the orientation of the multipole moments during the trajectory, conformational dependence of the atomic moments and stability of the simulations which are discussed here.

Characterization of molecular structures and properties of polyurethanes using molecular dynamics simulations

Thermotropic polyurethanes with mesogenic groups in side chains were prepared from two diisocyanates and four diols with stoichiometric ratios of reactive isocyanate (NCO) and hydroxy (OH) groups. Their thermal behavior was determined by differential scanning calorimetry. The effect of structure modifications of the diisocyanates and diols, in particular changes in the mesogen, were investigated. Introduction of mesogenic segments into the polymers suppresses the ordering. Stiff end substituents (phenyl and alkoxy groups) of the mesogens stabilize the mesophases to such an extent that the negative influence of long polymer chains is compensated and the liquid-crystalline properties are recovered. All-atom molecular dynamics simulations in the Cerius² modeling environment were carried out to characterize the structures of the polymers. Analysis of the dynamic trajectories at 20, 100, 120 and 170 °C revealed changes in conformation of macromolecules, which correlate with DSC measurements.Figure Example of structure relaxation of D4/TDI molecule at indicated simulation times (temperature 20 °C): a complete structure; b backbone structure; c top view of molecule     

Temperature dependence of protein-hydration hydrodynamics by molecular dynamics simulations

The dynamics of water molecules near the protein surface are different from those of bulk water and influence the structure and dynamics of the protein itself. To elucidate the temperature dependence hydration dynamics of water molecules, we present results from the molecular dynamic simulation of the water molecules surrounding two proteins (Carboxypeptidase inhibitor and Ovomucoid) at seven different temperatures (T=273 to 303 K, in increments of 5 K). Translational diffusion coefficients of the surface water and bulk water molecules were estimated from 2 ns molecular dynamics simulation trajectories. Temperature dependence of the estimated bulk water diffusion closely reflects the experimental values, while hydration water diffusion is retarded significantly due to the protein. Protein surface induced scaling of translational dynamics of the hydration waters is uniform over the temperature range studied, suggesting the importance protein-water interactions.     

Computer simulations of protein folding by targeted molecular dynamics

We have performed 128 folding and 45 unfolding molecular dynamics runs of chymotrypsin inhibitor 2 (CI2) with an implicit solvation model for a total simulation time of 0.4 microseconds. Folding requires that the three-dimensional structure of the native state is known. It was simulated at 300 K by supplementing the force field with a harmonic restraint which acts on the root-mean-square deviation and allows to decrease the distance to the target conformation. High temperature and/or the harmonic restraint were used to induce unfolding. Of the 62 folding simulations started from random conformations, 31 reached the native structure, while the success rate was 83% for the 66 trajectories which began from conformations unfolded by high-temperature dynamics. A funnel-like energy landscape is observed for unfolding at 475 K, while the unfolding runs at 300 K and 375 K as well as most of the folding trajectories have an almost flat energy landscape for conformations with less than about 50% of native contacts formed. The sequence of events, i.e., secondary and tertiary structure formation, is similar in all folding and unfolding simulations, despite the diversity of the pathways. Previous unfolding simulations of CI2 performed with different force fields showed a similar sequence of events. These results suggest that the topology of the native state plays an important role in the folding process.     

Study of the electrostatics treatment in molecular dynamics simulations

This article considers the treatment of long-range interactions in molecular dynamics simulations. We investigate the effects of using different cutoff distances, constant versus distance-dependent dielectric, and different smoothing methods. In contrast to findings of earlier studies, we find that increasing the cutoff over 8 A does not significantly improve the accuracy (Arnold and Ornstein, Proteins 1994;18:19-33), and using a distance-dependent dielectric instead of a constant dielectric also does not improve accuracy (Guenot and Kollman, Protein Sci 1992;1:1185-1205). This might depend on differences in simulation protocols or force fields, or both, because we use the CHARMM22 force field with stochastic boundary conditions, whereas earlier studies used other protocols and energy functions. We also note that the stability of the simulations is highly dependent on the starting structure, showing that accurate molecular simulations not only depend on a realistic simulation protocol but also on correct initial conditions.     

Understanding beta-hairpin formation by molecular dynamics simulations of unfolding

We have studied the mechanism of formation of a 16-residue beta-hairpin from the protein GB1 using molecular dynamics simulations in an aqueous environment. The analysis of unfolding trajectories at high temperatures suggests a refolding pathway consisting of several transient intermediates. The changes in the interaction energies of residues are related with the structural changes during the unfolding of the hairpin. The electrostatic energies of the residues in the turn region are found to be responsible for the transition between the folded state and the hydrophobic core state. The van der Waals interaction energies of the residues in the hydrophobic core reflect the behavior of the radius of gyration of the core region. We have examined the opposing influences of the protein-protein (PP) energy, which favors the native state, and the protein-solvent (PS) energy, which favors unfolding, in the formation of the beta-hairpin structure. It is found that the behavior of the electrostatic components of PP and PS energies reflects the structural changes associated with the loss of backbone hydrogen bonding. Relative changes in the PP and PS van der Waals interactions are related with the disruption of the hydrophobic core of a protein. The results of the simulations support the hydrophobic collapse mechanism of beta-hairpin folding.