Effect of the reaction field electrostatic term on the molecular dynamics simulation of the activation domain of procarboxypeptidase B

Molecular dynamics simulations of the activation domain of porcine procarboxypeptidase B (ADBp) were performed in order to examine the effects of the inclusion of a reaction field (RF) term into the calculation of electrostatics forces for highly charged proteins. Two simulations were performed with the GROMOS96 package, studying the influence of counterions on the final results. Comparison with previous results without the inclusion of the RF term (Martí-Renom, M.A., Mas,J.M., Oliva,B., Querol,E. and Avilés,F.X., Protein Engng, 1998, 11, 101-110) shows that the structure is well maintained when the RF term is included. Moreover, the analysis of the trajectories shows that simulations of solvated highly-charged proteins are sensitive to the presence of counterions, the secondary structures being more stable when their charges are neutralized.     

Importance of explicit salt ions for protein stability in molecular dynamics simulation.

The accurate and efficient treatment of electrostatic interactions is one of the challenging problems of molecular dynamics simulation. Truncation procedures such as switching or shifting energies or forces lead to artifacts and significantly reduced accuracy. The particle mesh Ewald (PME) method is one approach to overcome these problems by providing a computationally efficient means of calculating all long-range electrostatic interactions in a periodic simulation box by use of fast Fourier transformation techniques. For the application of the PME method to the simulation of a protein with a net charge in aqueous solution, counterions are added to neutralize the system. The usual procedure is to add charge-balancing counterions close to charged residues to neutralize the protein surface. In the present article, we show that for MD simulation of a small protein of marginal stability, the YAP-WW domain, explicit modeling of 0.2 M ionic strength (in addition to the charge-balancing counterions) is necessary to maintain a stable protein structure. Without explicit ions throughout the periodic simulation box, the charge-balancing counterions on the protein surface diffuse away from the protein, resulting in destruction of the beta-sheet secondary structure of the WW domain.     

Impact of Cl- and Na+ ions on simulated structure and dynamics of betaARK1 PH domain

A nonzero net charge of proteins at pH 7 is usually compensated by the addition of charge-balancing counter ions during molecular dynamics simulation, which reduces electrostatic interactions. For highly charged proteins, like the betaARK1 PH domain used here, it seems reasonable to also add explicit salt ions. To assess the impact of explicit salt ions, two molecular dynamics simulations of solvated betaARK1 PH domain have been carried out with different numbers of Cl- and Na+ ions, based on the Cornell et al. force field and the Ewald summation, which was used in the treatment of long-range electrostatic interactions. Initial positions of ions were obtained from the AMBER CION program. Increasing the number of ions alters the average structure in loop regions, as well as the fluctuation amplitudes of dihedral angles. We found unnaturally strong interactions between side chains in the absence of salt ions. The presence of salt ions reduces these electrostatic interactions. The time needed for the equilibration of the ionic environment around the protein, after initial placement of ions close to oppositely charged side chains, is in the nanosecond time range, which can be shortened by using a higher ionic strength. Our results also suggest selecting those methods that do not place the ions initially close to the protein surface.     

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.     

Molecular dynamics simulation of the unfolding of the human prion protein domain under low pH and high temperature conditions

Four 10-ns molecular dynamics (MD) simulations of the human prion protein domain (HuPrP 125-228) in explicit water solution have been performed. Each of the simulations mimicked a different environment of the protein: the neutral pH environment was simulated with all histidine residues neutral and bearing a ND proton and with other titratable side chains charged, the weakly acidic environment was simulated with all titratable side chains charged, the strongly acidic environment was simulated with all titratable side chains protonated. The protein in neutral pH environment was simulated at both ambient (298 K) and higher (350 K) temperatures. The native fold is stable in the neutral pH/ambient temperature simulation. Through out all other simulations, a quite stable core consisted of 10-20 residues around the disulfide bond retain their initial conformations. However, the secondary structures of the protein show changes of various degrees compared to the native fold, parts of the helices unfolded and the beta-sheets extended. Our simulations indicated that the heat-induced unfolding and acid-induced unfolding of HuPrP might follow different pathways: the initial stage of the acid-induced unfolding may include not only changes in secondary structures, but also changes in the tertiary structures. Under the strongly acidic condition, obvious tertiary structure changes take place after 10-ns simulation, the secondary structure elements and the loops becoming more parallel to each other, resulting in a compact state, which was stabilized by a large number of new, non-native side chain-side chain contacts. Such tertiary structure changes were not observed in the higher temperature simulation, and intuitively, they may favor the further extension of the beta-sheets and eventually the agglomeration of multiple protein molecules. The driving forces for this tertiary structure changes are discussed. Two additional 10-ns MD simulations, one with Asp202 protonated and the other with Glu196 protonated compared to the neutral pH simulation, were carried out. The results showed that the stability of the native fold is very subtle and can be strongly disturbed by eliminating a single negative charge at one of such key sites. Correlations of our results with previous experimental and theoretical studies are discussed.     

Prediction of titration properties of structures of a protein derived from molecular dynamics trajectories.

This paper explores the dependence of the molecular dynamics (MD) trajectory of a protein molecule on the titration state assigned to the molecule. Four 100-ps MD trajectories of bovine pancreatic trypsin inhibitor (BPTI) were generated, starting from two different structures, each of which was held in two different charge states. The two starting structures were the X-ray crystal structure and one of the solution structures determined by NMR, and the charge states differed only in the ionization state of N terminus. Although it is evident that the MD simulations were too short to sample fully the equilibrium distribution of structures in each case, standard Poisson-Boltzmann titration state analysis of the resulting configurations shows general agreement between the overall titration behavior of the protein and the charge state assumed during MD simulation: at pH 7, the total net charge of the protein resulting from the titration analysis is consistently lower for the protein with the N terminus assumed to be neutral than for the protein with the N terminus assumed to be charged. For most of the ionizable residues, the differences in the calculated pKaS among the four trajectories are statistically negligible and remain in good agreement with the data obtained by crystal structure titration and by experiment. The exceptions include the N terminus, which responds directly to the change of its imposed charge; the C terminus, which in the NMR structure interacts strongly with the former; and a few other residues (Arg 1, Glu 7, Tyr 35, and Arg 42) whose pKaS reflect the initial structure and the limited trajectory lengths. This study illustrates the importance of the careful assignment of protonation states at the start of MD simulations and points to the need for simulation methods that allow for the variation of the protonation state in the calculation of equilibrium properties.     

Counterion charge density determines the position and plasticity of RNA folding transition states

The self-assembly of RNA structure depends on the interactions of counterions with the RNA and with each other. Comparison of various polyamines showed that the tertiary structure of the Tetrahymena ribozyme is more stable when the counterions are small and highly charged. By monitoring the folding kinetics of the ribozyme as a function of polyamine concentration, we now find that the charge density of the counterions determines the positions of the folding transition states. The transition state ensemble (TSE) between U and N moves away from the native state as the counterion valence and charge density increase, as predicted by the Hammond postulate. The TSE is broader and less structured when the RNA is refolded in polyamines rather than Mg2+. That the charge density of the counterions determines the plasticity of the TSE demonstrates the importance of interactions among condensed counterions for the self-assembly of RNA structures. We propose that the major barrier to RNA folding is dominated by entropy changes when counterion charge density is low and enthalpy differences when it is high.     

Implicit solvent models for flexible protein-protein docking by molecular dynamics simulation

The suitability of three implicit solvent models for flexible protein-protein docking by procedures using molecular dynamics simulation is investigated. The three models are (i) the generalized Born (GB) model implemented in the program AMBER6.0; (ii) a distance-dependent dielectric (DDD) model; and (iii) a surface area-dependent model that we have parameterized and call the NPSA model. This is a distance-dependent dielectric model modified by neutralizing the ionizable side-chains and adding a surface area-dependent solvation term. These solvent models were first tested in molecular dynamics simulations at 300 K of the native structures of barnase, barstar, segment B1 of protein G, and three WW domains. These protein structures display a range of secondary structure contents and stabilities. Then, to investigate the performance of the implicit solvent models in protein docking, molecular dynamics simulations of barnase/barstar complexation, as well as PIN1 WW domain/peptide complexation, were conducted, starting from separated unbound structures. The simulations show that the NPSA model has significant advantages over the DDD and GB models in maintaining the native structures of the proteins and providing more accurate docked complexes.     

Hydrogen bonding in helical polypeptides from molecular dynamics simulations and amide hydrogen exchange analysis: alamethicin and melittin in methanol.

Molecular dynamics simulations of ion channel peptides alamethicin and melittin, solvated in methanol at 27 degrees C, were run with either regular alpha-helical starting structures (alamethicin, 1 ns; melittin 500 ps either with or without chloride counterions), or with the x-ray crystal coordinates of alamethicin as a starting structure (1 ns). The hydrogen bond patterns and stabilities were characterized by analysis of the dynamics trajectories with specified hydrogen bond angle and distance criteria, and were compared with hydrogen bond patterns and stabilities previously determined from high-resolution NMR structural analysis and amide hydrogen exchange measurements in methanol. The two alamethicin simulations rapidly converged to a persistent hydrogen bond pattern with a high level of 3(10) hydrogen bonding involving the amide NH's of residues 3, 4, 9, 15, and 18. The 3(10) hydrogen bonds stabilizing amide NH's of residues C-terminal to P2 and P14 were previously proposed to explain their high amide exchange stabilities. The absence, or low levels of 3(10) hydrogen bonds at the N-terminus or for A15 NH, respectively, in the melittin simulations, is also consistent with interpretations from amide exchange analysis. Perturbation of helical hydrogen bonding in the residues before P14 (Aib10-P14, alamethicin; T11-P14, melittin) was characterized in both peptides by variable hydrogen bond patterns that included pi and gamma hydrogen bonds. The general agreement in hydrogen bond patterns determined in the simulations and from spectroscopic analysis indicates that with suitable conditions (including solvent composition and counterions where required), local hydrogen-bonded secondary structure in helical peptides may be predicted from dynamics simulations from alpha-helical starting structures. Each peptide, particularly alamethicin, underwent some large amplitude structural fluctuations in which several hydrogen bonds were cooperatively broken. The recovery of the persistent hydrogen bonding patterns after these fluctuations demonstrates the stability of intramolecular hydrogen-bonded secondary structure in methanol (consistent with spectroscopic observations), and is promising for simulations on extended timescales to characterize the nature of the backbone fluctuations that underlie amide exchange from isolated helical polypeptides.     

Homology modelling and molecular dynamics simulations: comparative studies of human aquaporin-1

The structures of the mammalian water transport protein Aqp1 and of its bacterial homologue GlpF enables us to test whether homology models can be used to explore relationships between structure, dynamics and function in mammalian transport proteins. Molecular dynamics simulations (totalling almost 40 ns) were performed starting from: the X-ray structure of Aqp1; a homology model of Aqp1 based on the GlpF structure; and intermediate resolution structures of Aqp1 derived from electron microscopy. Comparisons of protein RMSDs vs. time suggest that the homology models are of comparable conformational stability to the X-ray structure, whereas the intermediate resolution structures exhibit significant conformation drift. For simulations based on the X-ray structure and on homology models, the flexibility profile vs. residue number correlates well with the crystallographic B-values for each residue. In the simulations based on intermediate resolution structures, mobility of the highly conserved NPA loops is substantially higher than in the simulations based on the X-ray structure or the homology models. Pore radius profiles remained relatively constant in the X-ray and homology model simulations but showed substantial fluctuations (reflecting the higher NPA loop mobility) in the intermediate resolution simulations. The orientation of the dipoles of water molecules within the pore is of key importance in maintaining low proton permeability through Aqp1. This property seems to be quite robust to the starting model used in the simulation. These simulations suggest that homology models based on bacterial homologues may be used to derive functionally relevant information on the structural dynamics of mammalian transport proteins.     

Protein structure and dynamics in nonaqueous solvents: insights from molecular dynamics simulation studies

Protein structure and dynamics in nonaqueous solvents are here investigated using molecular dynamics simulation studies, by considering two model proteins (ubiquitin and cutinase) in hexane, under varying hydration conditions. Ionization of the protein groups is treated assuming "pH memory," i.e., using the ionization states characteristic of aqueous solution. Neutralization of charged groups by counterions is done by considering a counterion for each charged group that cannot be made neutral by establishing a salt bridge with another charged group; this treatment is more physically reasonable for the nonaqueous situation, contrasting with the usual procedures. Our studies show that hydration has a profound effect on protein stability and flexibility in nonaqueous solvents. The structure becomes more nativelike with increasing values of hydration, up to a certain point, when further increases render it unstable and unfolding starts to occur. There is an optimal amount of water, approximately 10% (w/w), where the protein structure and flexibility are closer to the ones found in aqueous solution. This behavior can explain the experimentally known bell-shaped dependence of enzyme catalysis on hydration, and the molecular reasons for it are examined here. Water and counterions play a fundamental and dynamic role on protein stabilization, but they also seem to be important for protein unfolding at high percentages of bound water.     

Molecular dynamics simulation of human interleukin-4: comparison with NMR data and effect of pH,counterions and force field on tertiary structure stability

The human protein interleukin-4 (IL-4) has been simulated at two different pH values 2 and 6, with different amounts of counterions present in the aqueous solution, and with two different force-field parameter sets using molecular dynamics simulation with the aim of validation of force field and simulation set-up by comparison to experimental nuclear magnetic resonance data, such as proton–proton nuclear Overhauser effect (NOE) distance bounds, ³ J(HN,HC_α) coupling constants and backbone N–H order parameters. Thirteen simulations varying in the length from 3 to 7 ns are compared.

At pH 6 both force-field parameter sets used do largely reproduce the NOE's and order parameters, the GROMOS 45A3 set slightly better than the GROMOS 53A6 set. ³ J values predicted from the simulation agree less well with experimental values. At pH 2 the protein unfolds, unless counterions are explicitly present in the system, but even then the agreement with experiment is worse than at pH 6. When simulating a highly charged protein, such as IL-4 at pH 2, the inclusion of counterions in the simulation seems mandatory.     

Modeling H3 histone N-terminal tail and linker DNA interactions

Molecular dynamics computer simulations were performed for the 25-residue N-terminal tail of the H3 histone protein in the proximity of a DNA segment of 10 base pairs (bp), representing a model for the linker DNA in chromatin. Several least biased configurations were used as initial configurations. The secondary structure content of the protein was increased by the presence of DNA close to it, but the locations of the secondary motifs were different for different initial orientations of the DNA grooves with respect to the protein. As a common feature to all simulations, the electrostatic attraction between negatively charged DNA and positively charged protein was screened by the water solvent and counterbalanced by the intrinsic compaction of the protein due to hydrophobic effects. The protein secondary structure limited the covering of DNA by the protein to 4-5 bp. The degree of compaction and charge density of the bound protein suggests a possible role of H3 tail in a nonspecific bending and plasticity of the linker DNA when the protein is located in the crowded dense chromatin.     

Dynameomics: large-scale assessment of native protein flexibility

Structure is only the first step in understanding the interactions and functions of proteins. In this paper, we explore the flexibility of proteins across a broad database of over 250 solvated protein molecular dynamics simulations in water for an aggregate simulation time of approximately 6 micros. These simulations are from our Dynameomics project, and these proteins represent approximately 75% of all known protein structures. We employ principal component analysis of the atomic coordinates over time to determine the primary axis and magnitude of the flexibility of each atom in a simulation. This technique gives us both a database of flexibility for many protein fold families and a compact visual representation of a particular protein's native-state conformational space, neither of which are available using experimental methods alone. These tools allow us to better understand the nature of protein motion and to describe its relationship to other structural and dynamical characteristics. In addition to reporting general properties of protein flexibility and detailing many dynamic motifs, we characterize the relationship between protein native-state flexibility and early events in thermal unfolding and show that flexibility predicts how a protein will begin to unfold. We provide evidence that fold families have conserved flexibility patterns, and family members who deviate from the conserved patterns have very low sequence identity. Finally, we examine novel aspects of highly inflexible loops that are as important to structural integrity as conventional secondary structure. These loops, which are difficult if not impossible to locate without dynamic data, may constitute new structural motifs.     

Deducing hydration sites of a protein from molecular dynamics simulations

Invariant water molecules that are of structural or functional importance to proteins are detected from their presence in the same location in different crystal structures of the same protein or closely related proteins. In this study we have investigated the location of invariant water molecules from MD simulations of ribonuclease A, HIV1-protease and Hen egg white lysozyme. Snapshots of MD trajectories represent the structure of a dynamic protein molecule in a solvated environment as opposed to the static picture provided by crystallography. The MD results are compared to an analysis on crystal structures. A good correlation is observed between the two methods with more than half the hydration sites identified as invariant from crystal structures featuring as invariant in the MD simulations which include most of the functionally or structurally important residues. It is also seen that the propensities of occupying the various hydration sites on a protein for structures obtained from MD and crystallographic studies are different. In general MD simulations can be used to predict invariant hydration sites when there is a paucity of crystallographic data or to complement crystallographic results.     

Explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide duplex: lattice-sum versus reaction-field electrostatics

Five long-timescale (10 ns) explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide dimer are performed using the GROMOS 45A4 force field and the simple-point-charge water model, in order to investigate the effect of the treatment of long-range electrostatic interactions as well as of the box shape and size on the structure and dynamics of the molecule (starting from an idealised B-DNA conformation). Long-range electrostatic interactions are handled using either a lattice-sum (LS) method (particle–particle–particle–mesh; one simulation performed within a cubic box) or a cutoff-based reaction-field (RF) method (four simulations, with long-range cutoff distances of 1.4 or 2.0 nm and performed within cubic or truncated octahedral periodic boxes). The overall double-helical structure, including Watson–Crick (WC) base-pairing, is well conserved in the simulation employing the LS scheme. In contrast, the WC base-pairing is nearly completely disrupted in the four simulations employing the RF scheme. These four simulations result in highly distorted compact (cutoff distance of 1.4 nm) or extended (cutoff distance of 2 nm) structures, irrespective of the shape and size of the computational box. These differences observed between the two schemes seem correlated with large differences in the radial distribution function between charged entities (backbone phosphate groups and sodium counterions) within the system.     

Dynamics of water and ions around DNA: What is so special about them?

Water around biomolecules is special for behaving strangely – both in terms of structure and dynamics, while ions are found to control various interactions in biomolecules such as DNA, proteins and lipids. The questions that how water and ions around these biomolecules behave in terms of their structure and dynamics, and how they affect the biomolecular functions have triggered tremendous research activities worldwide. Such activities not only unfolded important static and dynamic properties of water and ions around these biomolecules, but also provoked heated debate regarding their explanation and role in biological functions. DNA, being negatively charged, interacts strongly with the surrounding dipolar water and positively charged counterions, leading to complex electrostatic coupling of water and ions with the DNA. Recent time-resolved fluorescence Stokes shift experiments and related computer simulation studies from our and other laboratories have unfolded some unique dynamic characteristics of water and ions near different structures of DNA. These results are discussed here to showcase the specialty of water and ion dynamics around DNA.     

Setting up and running molecular dynamics simulations of membrane proteins

Molecular dynamics simulations have become a popular and powerful technique to study lipids and membrane proteins. We present some general questions and issues that should be considered prior to embarking on molecular dynamics simulation studies of membrane proteins and review common simulation methods. We suggest a practical approach to setting up and running simulations of membrane proteins, and introduce two new (related) methods to embed a protein in a lipid bilayer. Both methods rely on placing lipids and the protein(s) on a widely spaced grid and then 'shrinking' the grid until the bilayer with the protein has the desired density, with lipids neatly packed around the protein. When starting from a grid based on a single lipid structure, or several potentially different lipid structures (method 1), the bilayer will start well-packed but requires more equilibration. When starting from a pre-equilibrated bilayer, either pure or mixed, most of the structure of the bilayer stays intact, reducing equilibration time (method 2). The main advantages of these methods are that they minimize equilibration time and can be almost completely automated, nearly eliminating one time consuming step in MD simulations of membrane proteins.     

Molecular dynamic simulations of cisplatin- and oxaliplatin-d(GG) intrastand cross-links reveal differences in their conformational dynamics

Mismatch repair proteins, DNA damage-recognition proteins and translesion DNA polymerases discriminate between Pt-GG adducts containing cis-diammine ligands (formed by cisplatin (CP) and carboplatin) and trans-RR-diaminocyclohexane ligands (formed by oxaliplatin (OX)) and this discrimination is thought to be important in determining differences in the efficacy, toxicity and mutagenicity of these platinum anticancer agents. We have postulated that these proteins recognize differences in conformation and/or conformational dynamics of the DNA containing the adducts. We have previously determined the NMR solution structure of OX-DNA, CP-DNA and undamaged duplex DNA in the 5'-d(CCTCAGGCCTCC)-3' sequence context and have shown the existence of several conformational differences in the vicinity of the Pt-GG adduct. Here we have used molecular dynamics simulations to explore differences in the conformational dynamics between OX-DNA, CP-DNA and undamaged DNA in the same sequence context. Twenty-five 10 ns unrestrained fully solvated molecular dynamics simulations were performed starting from two different DNA conformations using AMBER v8.0. All 25 simulations reached equilibrium within 4 ns, were independent of the starting structure and were in close agreement with previous crystal and NMR structures. Our data show that the cis-diammine (CP) ligand preferentially forms hydrogen bonds on the 5' side of the Pt-GG adduct, while the trans-RR-diaminocyclohexane (OX) ligand preferentially forms hydrogen bonds on the 3' side of the adduct. In addition, our data show that these differences in hydrogen bond formation are strongly correlated with differences in conformational dynamics, specifically the fraction of time spent in different DNA conformations in the vicinity of the adduct, for CP- and OX-DNA adducts. We postulate that differential recognition of CP- and OX-GG adducts by mismatch repair proteins, DNA damage-recognition proteins and DNA polymerases may be due, in part, to differences in the fraction of time that the adducts spend in a conformation favorable for protein binding.