Long time dynamics of Met-enkephalin: comparison of explicit and implicit solvent models

Met-enkephalin is one of the smallest opiate peptides. Yet, its dynamical structure and receptor docking mechanism are still not well understood. The conformational dynamics of this neuron peptide in liquid water are studied here by using all-atom molecular dynamics (MD) and implicit water Langevin dynamics (LD) simulations with AMBER potential functions and the three-site transferable intermolecular potential (TIP3P) model for water. To achieve the same simulation length in physical time, the full MD simulations require 200 times as much CPU time as the implicit water LD simulations. The solvent hydrophobicity and dielectric behavior are treated in the implicit solvent LD simulations by using a macroscopic solvation potential, a single dielectric constant, and atomic friction coefficients computed using the accessible surface area method with the TIP3P model water viscosity as determined here from MD simulations for pure TIP3P water. Both the local and the global dynamics obtained from the implicit solvent LD simulations agree very well with those from the explicit solvent MD simulations. The simulations provide insights into the conformational restrictions that are associated with the bioactivity of the opiate peptide dermorphin for the delta-receptor.     

A fluid salt-bridging cluster and the stabilization of p53

p53 is a homotetrameric tumor suppressor protein that is found to be mutated in most human cancers. Some of these mutations, particularly mutations to R337, fall in the tetramerization domain and cause defects in tetramer formation leading to loss of function. Mutation to His at this site has been found to destabilize the tetramer in a pH-dependent fashion. In structures of the tetramerization domain determined by crystallography, R337 from one monomer makes a salt bridge with D352 from another monomer, apparently helping to stabilize the tetramer. Here we present molecular dynamics simulations of wild-type p53 and the R337His mutant at several different pH and salt conditions. We find that the 337-352 salt bridge is joined by two other charged side chains, R333 and E349. These four residues do not settle into a fixed pattern of salt bridging, but continue to exchange salt-bridging partners on the nanosecond time scale throughout the simulation. This unusual system of fluid salt bridging may explain the previous finding from alanine scanning experiments that R333 contributes significantly to protein stability, even though in the crystal structure it is extended outward into solvent. This extended conformation of R333 appears to be the result of a specific crystal contact and, this contact being absent in the simulation, R333 turns inward to join its interaction partners. When R337 is mutated to His but remains positively charged, it maintains the original interaction with D352, but the newly observed interaction with E349 is weakened, accounting for the reduced stability of R337H even under mildly acidic conditions. When this His is deprotonated, the interaction with D352 is also lost, accounting for the further destabilization observed under mildly alkaline conditions. Simulations were carried out using both explicit and implicit solvent models, and both displayed similar behavior of the fluid salt-bridging cluster, suggesting that implicit solvent models can capture at least the qualitative features of this phenomenon as well as explicit solvent. Simulations under strongly acidic conditions in implicit solvent displayed the beginnings of the unfolding process, a destabilization of the hydrophobic dimer-dimer interface. Computational alanine scanning using the molecular mechanics Poisson-Boltzmann surface area method showed significant correlation to experimental unfolding data for charged and polar residues, but much weaker correlation for hydrophobic residues.     

Conformational free energy maps for globobiose (alpha-D-Galp-(1-->4)-beta-D-Galp) in implicit and explicit aqueous solution

Four Ramachandran maps of the conformational potential of mean force (PMF) for the galactose disaccharide globobiose (alpha-D-Galp-(1-->4)-beta-D-Galp) were calculated in vacuum, explicit water, with a simple high dielectric constant and a distance-dependent dielectric coefficient, respectively. This simple model of the galactan alpha-(1-->4)-linkage is shown to be conformationally restricted, with only a small range of syn-phi/syn-psi conformations predominating at standard temperature and pressure. This has implications for the preferred conformation and chain dynamics of alpha-galactosides. In addition, comparison of the relevant PMF surfaces reveals the substitution of a high dielectric constant for explicit water solution to be a valid approximation for reproducing the minimum energy conformation of this glycosidic linkage.     

Using molecular dynamics simulations on crambin to evaluate the suitability of different continuum dielectric and hydrogen atom models for protein simulations

Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 100-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of 2 and the all-hydrogen model. Continuum dielectric models of 1, 1r, and 2r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreement than the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.     

Molecular Dynamics Simulations of a DMPC Bilayer Using Nonadditive Interaction Models

We present a polarizable force field based on the charge-equilibration formalism for molecular dynamics simulations of phospholipid bilayers. We discuss refinement of headgroup dihedral potential parameters to reproduce ab initio conformational energies of dimethylphosphate calculated at the MP2/cc-pVTZ level of theory. We also address the refinement of electrostatic and Lennard-Jones (van der Waals) parameters to reproduce ab initio polarizabilities and water interaction energies of dimethylphosphate and tetramethylammonium. We present results of molecular dynamics simulations of a solvated dimyristoylphosphatidylcholine bilayer using this polarizable force field as well as the nonpolarizable, fixed-charge CHARMM27 and CHARMM27r force fields for comparison. Calculated atomic and electron-density profiles, deuterium order parameters, and headgroup orientations are found to be consistent with previous simulations and with experiment. Polarizable interaction models for solvent and lipid exhibit greater water penetration into the lipid interior; this is due to the variation of water molecular dipole moment from a bulk value of 2.6 Debye to a value of 1.9 Debye in the membrane interior. The reduction in the electrostatic component of the desolvation free-energy penalty allows for greater water density. The surface dipole potential predicted by the polarizable model is 0.95 V compared to the value of 0.8 V based on nonpolarizable force-field calculations. Effects of inclusion of explicit polarization are discussed in relation to water dipole moment and varying charge distributions. Dielectric permittivity profiles for polarizable and nonpolarizable interactions exhibit subtle differences arising from the nature of the individual component parameterizations; for the polarizable force field, we obtain a bulk dielectric permittivity of 79, whereas the nonpolarizable force field plateaus at 97 (the value for pure TIP3P water). In the membrane interior, both models predict unit permittivities, with the polarizable models contributing from one to two more units due to the optical dielectric (high-frequency dipole fluctuations). This contribution is a step toward the continuing development of a CHARMM (Chemistry at Harvard Molecular Mechanics) polarizable force field for simulations of biomacromolecular systems.     

Molecular dynamics simulations of small peptides: dependence on dielectric model and pH

There has been much interest recently in the structure of small peptides in solution. A recent study by Bradley and co-workers [(1989) in Techniques of Protein Chemistry, Hugli, T.E., Ed., Academic Press, Orlando, FL, pp. 531-546; (1990) Journal of Molecular Biology, 215, pp. 607-622] describes a 17-residue peptide that is stable as a monomeric helix in aqueous solution at low pH, as determined by two-dimensional nmr and CD spectroscopy. They also have determined the helix content of the peptide as a function of pH using CD. We performed molecular dynamics simulations, with an empirical force field, of this peptide at low pH, with three different dielectric models: a linear distance-dependent dielectric function (epsilon = R); a modified form [J. Ramstein and R. Lavery (1988) Proceedings of the National Academy of Science, USA, Vol. 85, pp. 7231-7235] of the sigmoidal distance-dependent dielectric function of Hingerty and co-workers [(1985) Biopolymers, Vol. 24, pp. 427-439]; and epsilon = 1 with the peptide immersed in a bath of water molecules. We found that simulations with the sigmoidal dielectric function and the model with explicit water molecules resulted in average distances for particular interactions that were consistent with the experimental nmr results, with the sigmoidal function best representing the data. However, these models exhibited very different helix-stabilizing interactions. We also performed simulations using the sigmoidal function at moderate and high pH to compare to experimental determinations of the pH dependence of helix content. Helix content did not decrease with increases in pH, as shown experimentally. We did, however, observe changes in a specific side chain-helix dipole interaction that was implicated in determining the pH-dependent behavior of this peptide. Overall, the sigmoidal dielectric function was a reasonable alternative to adding explicit water molecules. In comparing 100 ps molecular dynamics simulations, the sigmoidal function was much less computer intensive and sampled more of conformational space than the treatment using explicit water molecules. Sampling is especially important for this system since the peptide has been shown experimentally to populate both helical and nonhelical conformations.     

Molecular dynamics studies of a DNA-binding protein: 2. An evaluation of implicit and explicit solvent models for the molecular dynamics simulation of the Escherichia coli trp repressor.

Although aqueous simulations with periodic boundary conditions more accurately describe protein dynamics than in vacuo simulations, these are computationally intensive for most proteins. Trp repressor dynamic simulations with a small water shell surrounding the starting model yield protein trajectories that are markedly improved over gas phase, yet computationally efficient. Explicit water in molecular dynamics simulations maintains surface exposure of protein hydrophilic atoms and burial of hydrophobic atoms by opposing the otherwise asymmetric protein-protein forces. This properly orients protein surface side chains, reduces protein fluctuations, and lowers the overall root mean square deviation from the crystal structure. For simulations with crystallographic waters only, a linear or sigmoidal distance-dependent dielectric yields a much better trajectory than does a constant dielectric model. As more water is added to the starting model, the differences between using distance-dependent and constant dielectric models becomes smaller, although the linear distance-dependent dielectric yields an average structure closer to the crystal structure than does a constant dielectric model. Multiplicative constants greater than one, for the linear distance-dependent dielectric simulations, produced trajectories that are progressively worse in describing trp repressor dynamics. Simulations of bovine pancreatic trypsin were used to ensure that the trp repressor results were not protein dependent and to explore the effect of the nonbonded cutoff on the distance-dependent and constant dielectric simulation models. The nonbonded cutoff markedly affected the constant but not distance-dependent dielectric bovine pancreatic trypsin inhibitor simulations. As with trp repressor, the distance-dependent dielectric model with a shell of water surrounding the protein produced a trajectory in better agreement with the crystal structure than a constant dielectric model, and the physical properties of the trajectory average structure, both with and without a nonbonded cutoff, were comparable.     

Molecular dynamics effects on protein electrostatics

Electrostatic calculations have been carried out on a number of structural conformers of tuna cytochrome c. Conformers were generated using molecular dynamics simulations with a range of solvent simulating, macroscopic dielectric formalisms, and one solvent model that explicitly included solvent water molecules. Structures generated using the lowest dielectric models were relatively tight, with side chains collapsed on the surface, while those from the higher dielectric models had more internal and external fluidity, with surface side chains exploring a fuller range of conformational space. The average structure generated with the explicitly solvated model corresponded most closely with the crystal structure. Individual pK values, overall titration curves, and electrostatic potential surfaces were calculated for average structures and structures along each simulation. Differences between structural conformers within each simulation give rise to substantial changes in calculated local electrostatic interactions, resulting in pK value fluctuations for individual sites in the protein that vary by 0.3-2.0 pK units from the calculated time average. These variations are due to the thermal side chain reorientations that produce fluctuations in charge site separations. Properties like overall titration curves and pH dependent stability are not as sensitive to side chain fluctuations within a simulation, but there are substantial effects between simulations due to marked differences in average side chain behavior. These findings underscore the importance of proper dielectric formalism in molecular dynamics simulations when used to generate alternate solution structures from a crystal structure, and suggest that conformers significantly removed from the average structure have altered electrostatic properties that may prove important in episodic protein properties such as catalysis.     

Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations

Abstract

Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 100-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of 2 and the all-hydrogen model. Continuum dielectric models of 1,1^*r, and 2^*r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreement than the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4^*r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Molecular dynamics simulations of the cytochrome c3-rubredoxin complex from Desulfovibrio vulgaris.

Molecular dynamics simulations have been carried out on the complex formed between the tetraheme cytochrome c3 and the iron protein rubredoxin from the sulfate-reducing bacterium Desulfovibrio vulgaris. These simulations were performed both with explicit solvent water molecules included, and without solvent molecules using a distance-dependent dielectric constant to approximate the screening effects of solvent. The results of both simulations are strikingly different, indicating that the representation of environmental effects is important in such simulations. For example, a striking adaptation of the two proteins seen in the nonsolvated simulation is not seen when explicit solvent water is included; in fact, the complex appears to become weaker in the solvated simulation. Nonetheless, the iron-iron distance decreases more significantly in the solvated simulation than in the nonsolvated simulation. It was found that in both cases molecular dynamics optimized the structures further than energy minimization alone.     

Structure and dynamics of human vimentin intermediate filament dimer and tetramer in explicit and implicit solvent models

Intermediate filaments, in addition to microtubules and microfilaments, are one of the three major components of the cytoskeleton in eukaryotic cells, and play an important role in mechanotransduction as well as in providing mechanical stability to cells at large stretch. The molecular structures, mechanical and dynamical properties of the intermediate filament basic building blocks, the dimer and the tetramer, however, have remained elusive due to persistent experimental challenges owing to the large size and fibrillar geometry of this protein. We have recently reported an atomistic-level model of the human vimentin dimer and tetramer, obtained through a bottom-up approach based on structural optimization via molecular simulation based on an implicit solvent model (Qin et al. in PLoS ONE 2009 4(10):e7294, 9). Here we present extensive simulations and structural analyses of the model based on ultra large-scale atomistic-level simulations in an explicit solvent model, with system sizes exceeding 500,000 atoms and simulations carried out at 20 ns time-scales. We report a detailed comparison of the structural and dynamical behavior of this large biomolecular model with implicit and explicit solvent models. Our simulations confirm the stability of the molecular model and provide insight into the dynamical properties of the dimer and tetramer. Specifically, our simulations reveal a heterogeneous distribution of the bending stiffness along the molecular axis with the formation of rather soft and highly flexible hinge-like regions defined by non-alpha-helical linker domains. We report a comparison of Ramachandran maps and the solvent accessible surface area between implicit and explicit solvent models, and compute the persistence length of the dimer and tetramer structure of vimentin intermediate filaments for various subdomains of the protein. Our simulations provide detailed insight into the dynamical properties of the vimentin dimer and tetramer intermediate filament building blocks, which may guide the development of novel coarse-grained models of intermediate filaments, and could also help in understanding assembly mechanisms.     

Effective water model for Monte Carlo simulations of proteins

We present an effective theory for water. Our goal is to formulate on accurate model for the effects of solvation on protein dynamics, without incurring the huge computational cost and the slow temporal evolution typical of molecular dynamics simulations of liquids. We replace the individual water molecules in an all-atom potential with a local dielectric density field, with self interactions given by the Landau-Ginzburg free energy and external interactions by Lennard-Jones forces at the surface of the protein atoms. We explore conformational space with finite temperature Monte Carlo dynamics, using parallel Langevin and Fourier acceleration algorithms well suited to data-parallel computer architectures such as the Connection Machine. To establish the validity of our approximations, we compare our electrostatic contribution to the solvalion energy with the results of Lim, Bashford, and Karplus using a conventional static continuum dielectric cavity model, and the non electrostatic contributions with estimates of hydrophohic surface free energy. Our model can also accommodate ionic charges and temperature fluctuations, We propose future investigations extending our effective theory of solvation to include explicit orientational entropy and hydroxen-bonding terms. © 1995 John Wiley & Sons, Inc.     

Pathways of H2 toward the active site of [NiFe]-hydrogenase

Hydrogenases catalyze the reversible oxidation of molecular hydrogen (H(2)), but little is known about the diffusion of H(2) toward the active site. Here we analyze pathways for H(2) permeation using molecular dynamics (MD) simulations in explicit solvent. Various MD simulation replicates were done, to improve the sampling of the system states. H(2) easily permeates hydrogenase in every simulation and it moves preferentially in channels. All H(2) molecules that reach the active site made their approach from the side of the Ni ion. H(2) is able to reach distances of <4 A from the active site, although after 6 A permeation is difficult. In this region we mutated Val-67 into alanine and perform new MD simulations. These simulations show an increase of H(2) inside the protein and at lower distances from the active site. This valine can be a control point in the H(2) access to the active center.     

Unbinding free energy of acetylcholinesterase bound oxime drugs along the gorge pathway from metadynamics‐umbrella sampling investigation

Because of the pivotal role that the nerve enzyme, acetylcholinesterase plays in terminating nerve impulses at cholinergic synapses. Its active site, located deep inside a 20 Å gorge, is a vulnerable target of the lethal organophosphorus compounds. Potent reactivators of the intoxicated enzyme are nucleophiles, such as bispyridinium oxime that binds to the peripheral anionic site and the active site of the enzyme through suitable cation–π interactions. Atomic scale molecular dynamics and free energy calculations in explicit water are used to study unbinding pathways of two oxime drugs (Ortho‐7 and Obidoxime) from the gorge of the enzyme. The role of enzyme‐drug cation–π interactions are explored with the metadynamics simulation. The metadynamics discovered potential of mean force (PMF) of the unbinding events is refined by the umbrella sampling (US) corrections. The bidimensional free energy landscape of the metadynamics runs are further subjected to finite temperature string analysis to obtain the transition tube connecting the minima and bottlenecks of the unbinding pathway. The PMF is also obtained from US simulations using the biasing potential constructed from the transition tube and are found to be consistent with the metadynamics‐US corrected results. Although experimental structural data clearly shows analogous coordination of the two drugs inside the gorge in the bound state, the PMF of the drug trafficking along the gorge pathway point, within an equilibrium free energy context, to a multistep process that differs from one another. Routes, milestones and subtlety toward the unbinding pathway of the two oximes at finite temperature are identified. Proteins 2014; 82:1799–1818. © 2014 Wiley Periodicals, Inc.     

An evaluation of implicit and explicit solvent model systems for the molecular dynamics simulation of bacteriophage T4 lysozyme

In this report we examine several solvent models for use in molecular dynamics simulations of protein molecules with the Discover program from Biosym Technologies. Our goal was to find a solvent system which strikes a reasonable balance among theoretical rigor, computational efficiency, and experimental reality. We chose phage T4 lysozyme as our model protein and analyzed 14 simulations using different solvent models. We tested both implicit and explicit solvent models using either a linear distance-dependent dielectric or a constant dielectric. Use of a linear distance-dependent dielectric with implicit solvent significantly diminished atomic fluctuations in the protein and kept the protein close to the starting crystal structure. In systems using a constant dielectric and explicit solvent, atomic fluctuations were much greater and the protein was able to sample a larger portion of conformational space. A series of nonbonded cutoff distances (9.0, 11.5, 15.0, 20.0 Å) using both abrupt and smooth truncation of the nonbonded cutoff distances were tested. The method of dual cutoffs was also tested. We found that a minimum nonbonded cutoff distance of 15.0 Å was needed in order to properly couple solvent and solute. Distances shorter than 15.0 Å resulted in a significant temperature gradient between the solvent and solute. In all trajectories using the proprietary Discover switching function, we found significant denaturation in the protein backbone; we were able to run successful trajectories only in those simulations that used no switching function. We were able to significantly reduce the computational burden by using dual cutoffs and still calculate a quality trajectory. In this method, we found that an outer cutoff distance of 15.0 Å and an inner cutoff distance of 11.5 worked well. While a 10 Å shell of explicit water yielded the best results, a 6 A shell of water yielded satisfactory results with nearly a 40% reduction in computational cost. © 1994 John Wiley & Sons, Inc.     

Open-Boundary Molecular Dynamics of a DNA Molecule in a Hybrid Explicit/Implicit Salt Solution

The composition and electrolyte concentration of the aqueous bathing environment have important consequences for many biological processes and can profoundly affect the behavior of biomolecules. Nevertheless, because of computational limitations, many molecular simulations of biophysical systems can be performed only at specific ionic conditions: either at nominally zero salt concentration, i.e., including only counterions enforcing the system’s electroneutrality, or at excessive salt concentrations. Here, we introduce an efficient molecular dynamics simulation approach for an atomistic DNA molecule at realistic physiological ionic conditions. The simulations are performed by employing the open-boundary molecular dynamics method that allows for simulation of open systems that can exchange mass and linear momentum with the environment. In our open-boundary molecular dynamics approach, the computational burden is drastically alleviated by embedding the DNA molecule in a mixed explicit/implicit salt-bathing solution. In the explicit domain, the water molecules and ions are both overtly present in the system, whereas in the implicit water domain, only the ions are explicitly present and the water is described as a continuous dielectric medium. Water molecules are inserted and deleted into/from the system in the intermediate buffer domain that acts as a water reservoir to the explicit domain, with both water molecules and ions free to enter or leave the explicit domain. Our approach is general and allows for efficient molecular simulations of biomolecules solvated in bathing salt solutions at any ionic strength condition.     

Molecular dynamics simulations for the prediction of the dielectric spectra of alcohols,glycols and monoethanolamine

The response of molecular systems to electromagnetic radiation in the microwave region (0.3–300 GHz) has been principally studied experimentally, using broadband dielectric spectroscopy. However, relaxation times corresponding to reorganisation of molecular dipoles due to their interaction with electromagnetic radiation at microwave frequencies are within the scope of modern molecular simulations. In this work, fluctuations of the total dipole moment of a molecular system, obtained through molecular dynamics simulations, are used to determine the dielectric spectra of water, a series of alcohols and glycols, and monoethanolamine. Although the force fields employed in this study have principally been developed to describe thermodynamic properties, most them give fairly good predictions of this dynamical property for these systems. However, the inaccuracy of some models and the long simulation times required for the accurate estimation of the static dielectric constant can sometimes be problematic. We show that the use of the experimental value for the static dielectric constant in the calculations, instead of the one predicted by the different models, yields satisfactory results for the dielectric spectra, and hence the heat absorbed from microwaves, avoiding the need for extraordinarily long simulations or re-calibration of molecular models.     

Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

A continuum electrostatics approach for molecular dynamics (MD) simulations of macromolecules is presented and analyzed for its performance on a peptide and a globular protein. The approach incorporates the screened Coulomb potential (SCP) continuum model of electrostatics, which was reported earlier. The model was validated in a broad set of tests some of which were based on Monte Carlo simulations that included single amino acids, peptides, and proteins. The implementation for large-scale MD simulations presented in this article is based on a pairwise potential that makes the electrostatic model suitable for fast analytical calculation of forces. To assess the suitability of the approach, a preliminary validation is conducted, which consists of (i) a 3-ns MD simulation of the immunoglobulin-binding domain of streptococcal protein G, a 56-residue globular protein and (ii) a 3-ns simulation of Dynorphin, a biological peptide of 17 amino acids. In both cases, the results are compared with those obtained from MD simulations using explicit water (EW) molecules in an all-atom representation. The initial structure of Dynorphin was assumed to be an alpha-helix between residues 1 and 9 as suggested from NMR measurements in micelles. The results obtained in the MD simulations show that the helical structure collapses early in the simulation, a behavior observed in the EW simulation and consistent with spectroscopic data that suggest that the peptide may adopt mainly an extended conformation in water. The dynamics of protein G calculated with the SCP implicit solvent model (SCP-ISM) reveals a stable structure that conserves all the elements of secondary structure throughout the entire simulation time. The average structures calculated from the trajectories with the implicit and explicit solvent models had a cRMSD of 1.1 A, whereas each average structure had a cRMSD of about 0.8A with respect to the X-ray structure. The main conformational differences of the average structures with respect to the crystal structure occur in the loop involving residues 8-14. Despite the overall similarity of the simulated dynamics with EW and SCP models, fluctuations of side-chains are larger when the implicit solvent is used, especially in solvent exposed side-chains. The MD simulation of Dynorphin was extended to 40 ns to study its behavior in an aqueous environment. This long simulation showed that the peptide has a tendency to form an alpha-helical structure in water, but the stabilization free energy is too weak, resulting in frequent interconversions between random and helical conformations during the simulation time. The results reported here suggest that the SCP implicit solvent model is adequate to describe electrostatic effects in MD simulation of both peptides and proteins using the same set of parameters. It is suggested that the present approach could form the basis for the development of a reliable and general continuum approach for use in molecular biology, and directions are outlined for attaining this long-term goal.