GROMOS 53A6 Research Articles

A Comparative Study for Molecular Dynamics Simulations of Liquid Benzene.

The classical equilibrium and nonequilibrium molecular dynamics simulations for liquid benzene, the prototypical aromatic π-π interaction system, are performed using a variety of molecular force fields, OPT-FF, AMBER 03, general AMBER force field (GAFF), OPLS-AA, OPLS-CS, CHARMM27, GROMOS 53A5, and GROMOS 53A6. The simulated results of the molecular structure and thermodynamic properties of liquid benzene are compared with the experimental data available in the literature, accounting for the superiority of each force field in the descriptions of the π-π interaction system. The OPLS-AA force field is recommended to be the best one, which reproduces quite well the properties examined in this work, while the others fail in predicting either the local structure or the thermodynamic properties. Such distinct discrepancies for the above force fields are discussed within the scheme of the pairwise interaction construction of the standard force field, which will stimulate searching for a force field with generally good quality not only in terms of microstructure descriptions but also in the predictions of the thermodynamic properties of the liquids.

Lipid models for united-atom molecular dynamics simulations of protein

AbstractMolecular dynamics (MD) simulations of biomolecules are performed at various levels of atomic detail from all-atom models to coarse grained models. United-atom models are at an intermediate level, in which non-polar hydrogen atoms are subsumed into their adjacent carbon-atom reducing the number of particles for a DPPC lipid molecule by 60%. United-atom models of the lipids DPPC, DMPC, POPC and POPG in the GROMOS96 53a6 force field were developed, that reproduced experimental properties of lipid bilayers without assumption of a constant surface area or inclusion of surface pressure. In the absence of experimental data a surface per lipid area for the negatively charged POPG of 0.700±0.007 nm2 was obtained that is higher than 0.53 to 0.56 nm2 reported in previous simulation studies in the literature. It is argued that MD simulations of anionic bilayers may have underestimated the steric requirements of sodium counter ions entering the head group region in order to compensate the negative charge. The use of the DMPC GROMOS96 53a6 model was evaluated in an MD simulation of the ErbB2 transmembrane peptide and showed closer agreement with the NMR-derived structure than the GROMOS87 force field that was most commonly used for membrane protein simulations previously. Lipid topologies and coordinates are available via LipidBook (http://lipidbook.bioch.ox.ac.uk/).

Lipid models for united-atom molecular dynamics simulations of protein

Molecular dynamics (MD) simulations of biomolecules are performed at various levels of atomic detail from all-atom models to coarse grained models. United-atom models are at an intermediate level, in which non-polar hydrogen atoms are subsumed into their adjacent carbon-atom reducing the number of particles for a DPPC lipid molecule by 60%. United-atom models of the lipids DPPC, DMPC, POPC and POPG in the GROMOS96 53a6 force field were developed, that reproduced experimental properties of lipid bilayers without assumption of a constant surface area or inclusion of surface pressure. In the absence of experimental data a surface per lipid area for the negatively charged POPG of 0.700±0.007 nm2 was obtained that is higher than 0.53 to 0.56 nm2 reported in previous simulation studies in the literature. It is argued that MD simulations of anionic bilayers may have underestimated the steric requirements of sodium counter ions entering the head group region in order to compensate the negative charge. The use of the DMPC GROMOS96 53a6 model was evaluated in an MD simulation of the ErbB2 transmembrane peptide and showed closer agreement with the NMR-derived structure than the GROMOS87 force field that was most commonly used for membrane protein simulations previously. Lipid topologies and coordinates are available via LipidBook (http://lipidbook.bioch.ox.ac.uk/).

Inhibition of Aggregation of Amyloid Peptides by Beta-Sheet Breaker Peptides and Their Binding Affinity

The effects of beta-sheet breaker peptides KLVFF and LPFFD on the oligomerization of amyloid peptides were studied by all-atom simulations. It was found that LPFFD interferes the aggregation of Aβ(16-22) peptides to a greater extent than does KLVFF. Using the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method, we found that the former binds more strongly to Aβ(16-22). Therefore, by simulations, we have clarified the relationship between aggregation rates and binding affinity: the stronger the ligand binding, the slower the oligomerization process. The binding affinity of pentapeptides to full-length peptide Aβ(1-40) and its mature fibrils has been considered using the Autodock and MM-PBSA methods. The hydrophobic interaction between ligands and receptors plays a more important role for association than does hydrogen bonding. The influence of beta-sheet breaker peptides on the secondary structures of monomer Aβ(1-40) was studied in detail, and it turns out that, in their presence, the total beta-sheet content can be enhanced. However, the aggregation can be slowed because the beta-content is reduced in fibril-prone regions. Both pentapeptides strongly bind to monomer Aβ(1-40), as well as to mature fibrils, but KLVFF displays a lower binding affinity than LPFFD. Our findings are in accord with earlier experiments that both of these peptides can serve as prominent inhibitors. In addition, we predict that LPFFD inhibits/degrades the fibrillogenesis of full-length amyloid peptides better than KLVFF. This is probably related to a difference in their total hydrophobicities in that the higher the hydrophobicity, the lower the inhibitory capacity. The GROMOS96 43a1 force field with explicit water and the force field proposed by Morris et al. (Morris et al. J. Comput. Chem. 1998, 19, 1639 ) were employed for all-atom molecular dynamics simulations and Autodock experiments, respectively.

Modification of Martini force field for molecular dynamics simulation of hydrophobic charge induction chromatography of lysozyme

Modeling, especially the force field, is crucial for the accuracy of molecular dynamics (MD) simulations. In order for more accurate description of adsorption and desorption behaviors of lysozyme in hydrophobic charge induction chromatography (HCIC), the Martini coarse-grained (CG) force field has been modified based on the statistical analysis and comparison of an all-atom (AA) force field, GROMOS96 43A1, and the Martini force field. The parameters describing the protein–adsorbent interactions have been adjusted to avoid too strong and unrealistic adsorption of lysozyme on the agarose matrix and HCIC ligands. It is found that the adsorption and desorption behaviors monitored using the modified Martini force field and MD simulation are consistent with previous simulation results with 46-bead β-barrel model protein. Repeated adjustment of both protein position and orientation is necessary to generate enough contacts for a stable adsorption. After reducing the pH in the mobile phase, the lysozyme–ligand electrostatic repulsion leads to protein desorption. In the adsorption process, little conformational transition of lysozyme is observed due to its stable structure, which is in line with previous experimental observations. So, it is concluded that after appropriate modification, the Martini force field can be used to examine the HCIC process of lysozyme. The modification strategy has thus extended the applicability of the Martini force field to protein chromatography, and it is expected to facilitate studies of exploring the molecular details in adsorption chromatography of proteins.

Headgroup Mediated Water Insertion into the DPPC Bilayer: A Molecular Dynamics Study

Molecular dynamics simulation was performed on the 1,2-dipalmitoyl-sn-phosphocholine (DPPC) bilayer-water system using the GROMOS96 53a6 united atom force field. The transferability of force field was tested by reproducing the area per lipid within 3% accuracy from the experimental value. The simulation shows that water can penetrate much deeper inside the bilayer almost up to the starting point of the aliphatic chain. There is significant evidence from experiments that water goes deep in the DPPC bilayer, but it has not been reported from theoretical work. The mechanism of insertion of water deep inside the lipid bilayer is still not clear. In this report, for the first time, the mechanism of water insertion deep into the bilayer has been proposed. Water transport occurs by the headgroup and its first solvation shell. The trimethyl ammonium (NMe(3)) group (headgroup of DPPC) has two stable conformations at the bilayer-water interface, one outside the bilayer and another inside it. The NMe(3) group has a large clustering of water around it and takes the water molecules inside the bilayer with it during its entry into the bilayer. The water molecules penetrate into the bilayer with the help of the NMe(3) group present at the headgroup of DPPC and eventually form hydrogen bonds with carbonyl oxygen present deep inside the bilayer. Structural characteristics at the bilayer-water interface region are also reported.

Atomistic and Coarse-Grain Computer Simulations of the Assembled AcrAB-TolC Multidrug Efflux Pump

With known antibiotics losing their efficiency faster than new ones can be developed, a better understanding of the underlying molecular mechanisms is paramount. Multidrug resistance is often caused by an over-production of efflux transporters that expel drug molecules before they can affect their targets inside the bacterial cell.In Escherichia coli, AcrAB-TolC serves as the major multidrug efflux pump using proton-motive force over the inner membrane to extrude drugs out of the cell. While X-ray structures have been solved separately for the individual components, the best structural information on the assembled efflux pump is a docking structure based on biochemical cross-linking data [1]. To gain insight into the conformational dynamics of AcrB, AcrA & TolC in complex, we have embedded the docking structure in two phospholipid bilayers solvated in a 150 mM NaCl solution. Asymmetric AcrB was considered and each monomer simulated in a different protonation state as suggested in [2]. To assess the influence of AcrA, we set up a second model without the periplasmic adaptor protein.To study subunit interplay, complex stability, flexibility, protein-membrane interactions and the dynamics of transport pathways, we carried out multiple molecular dynamics simulations of Acr(A)B-TolC in atomistic and coarse-grained representations on a nanosecond and microsecond time scale. All simulations were carried out in GroMACS 4.0 using the GROMOS96 53a6 and MARTINI force fields.[1] Symmons et al. (2008)-PNAS_106(17):7173–8[2] Pos (2009)-BBA_1794(5):782–93

Simulated Substrate Binding to the Inner Membrane Translocase AcrB

Functioning as the engine of the AcrAB-TolC efflux pump, AcrB plays a key role in multidrug resistance in Escherichia coli. AcrB occurs as a trimer, where each monomer has a different conformation: loose (L), tight (T) or open (O) - representing different consecutive states in the reaction cycle. Here we report molecular dynamics simulations of the asymmetric 2GIF AcrB x-ray structure embedded in a phospholipid bilayer / 150 mM NaCl environment. To study the interaction with one of AcrB's simplest substrates, 25 hexane molecules were added to the system, with three hexanes placed directly in front of each monomer's porter domain. Using GroMACS 4.0.3 and the GROMOS96 53a6 force field we performed 5 independent MD runs, each 50ns long. Three runs were further extended by 150 ns. During one simulation, we observe a hexane entering the presumed drug transport channel of the (L) monomer. Binding occurs in a stepwise process during which the hexane moves towards the hydrophobic binding pocket inside the protein before reaching a final position after 25ns, 9.7 Å away from Phe-628. For the (T) monomer we observe hexane binding in a single step, reaching its final position 9.2 Å away from Phe-628. In none of the runs substrate binding takes place in the (O) monomer. However, an accumulation of hexane molecules in front of the closed porter domain could be observed here. In a second approach, we repeated the simulations now providing the antibiotic chloramphenicol as transport substrate instead. Correction et al.Biophysical JournalMay 04, 2011In Brief2963-Pos. Full-Text PDF Open Archive

Effects of pH and Temperature on the Structural and Thermodynamic Character of a-syn12 Peptide in Aqueous Solution

The structural and thermodynamic characters of α-syn12 peptide in aqueous solution at different pH and temperatures have been investigated through temperature replica exchange molecular dynamics (T-REMD) simulations with GROMOS 43A1 force field. The two independent T-REMD simulations were completed at pH = 7.0 and 10.0, respectively. Each replica was run for 300 ns. The structural and thermodynamic characters of α-syn12 peptide were studied based on the distributions of backbone dihedral angles, the free energy surface, and the stability of different type structure and the favorite conformations of the peptide. The results showed that the simulation at pH = 10.0 produced more sampling in α region than the simulation at pH = 7.0. The temperature changes from 283 K to 308 K result in negligible effects on the distributions of backbone dihedral angle. The β hairpin conformation with Turn9-6 and four hydrogen bonds (HB4-11, HB6-9, HB9-6 and HB11-4) is the lowest free energy state in the simulation at pH = 7.0. However, for the simulation at pH = 10.0, the lowest free energy state corresponds to a structure with Turn9-6 and two hydrogen bonds (HB6-10 and HB10-6) induced by an overly strong residue-residue interaction effect between lysine residues. For the simulation at pH = 7.0, the free energy change of the α-syn12 peptide from the unfolded state to the β hairpin state was in good agreement with the experiments and molecular dynamics simulation results for the other β-peptides, the β hairpin state of the α-syn12 peptide included the conformations that not only the Turn9-6 is formed, but also the terminus are closed together in space. However, the subtle balances between lysine-lysine interactions and lysine-solvent interaction are disrupted in the simulation at pH = 10.0, which induced the assembly of lysine residues, the β hairpin conformation is destabilized by the deprotonation of the Lys side chain. This study can help us to understand the conformation changes and the thermodynamic character of α-syn12 peptide at atomic level induced by changing pH and temperature, which is propitious to reveal the nosogenesis of Parkinson disease. In our knowledge, this is the first report to study the influence of pH and temperature on isolated α-syn12 peptide in water by T-REMD.

Molecular dynamics and experimental study of conformation change of poly(N-isopropylacrylamide) hydrogels in water

The conformation change of poly(N-isopropylacrylamide) (PNIPAAm) hydrogel as a function of temperature is studied both experimentally and by molecular simulation. The experiments show a Θ -temperature of approximately 305 K and a width of the transition region of approximately 20 K, independent of the amount of cross-linker. The conformation change of the hydrogel is caused by a conformation change of the single polymer chains in its backbone. The latter was studied by massively parallel molecular dynamics simulations with GROMACS 4. Different force fields for describing the hydrogel backbone (GROMOS-87, GROMOS-96 53A6 and OPLS-AA) as well as the water molecule (SPC/E and TIP4P) were investigated. In the simulations the mean radius of gyration of the polymer chains was monitored. Two force field combinations ( GROMOS- 96 53 A 6 + TIP 4 P and OPLS-AA + SPC/E ) yield the conformation change upon variation of the temperature. The results for OPLS-AA + SPC/E are in fair agreement with the experimental results.

Amiloride Docking to Acid-sensing Ion Channel-1

Amiloride is a small molecule diuretic, which has been used to dissect sodium transport pathways in many different systems. This drug is known to interact with the epithelial sodium channel and acid-sensing ion channel proteins, as well as sodium/hydrogen antiporters and sodium/calcium exchangers. The exact structural basis for these interactions has not been elucidated as crystal structures of these proteins have been challenging to obtain, though some involved residues and domains have been mapped. This work examines the interaction of amiloride with acid-sensing ion channel-1, a protein whose structure is available using computational and experimental techniques. Using molecular docking software, amiloride and related molecules were docked to model structures of homomeric human ASIC-1 to generate potential interaction sites and predict which analogs would be more or less potent than amiloride. The predictions made were experimentally tested using whole-cell patch clamp. Drugs previously classified as NCX or NHE inhibitors are shown to also inhibit hASIC-1. Potential docking sites were re-examined against experimental data to remove spurious interaction sites. The voltage sensitivity of inhibitors was also examined. Using the aggregated data from these computational and experimental experiments, putative interaction sites for amiloride and hASIC-1 have been defined. Future work will experimentally verify these interaction sites, but at present this should allow for virtual screening of drug libraries at these putative interaction sites.

Mapping the Proton Conduction Pathways in the Inner Membrane Multi-Drug Translocase AcrB

With increasing frequency the emergence of new drug resistances constitutes a serious problem in cancer therapy or the treatment of bacterial infections. A major mechanism bacteria employ against antibiotics is based on multidrug efflux pumps extruding the drugs before they can reach their targets. In Escherichia coli the major efflux system comprises the inner-membrane translocase AcrB, the membrane-fusion protein AcrA and the outer-membrane channel TolC. AcrB functions as the engine of this complex, using proton motive force to expel a wide variety of unrelated toxic compounds such as antibiotics, disinfectants or detergents. The molecular details of how proton conduction through AcrB is coupled to drug expulsion are not fully understood yet. To gain insight into the AcrB proton conduction pathway we performed a series of 5 independent molecular dynamics simulations of AcrB in a phospholipid/water environment at a 150mMol NaCl concentration. Each monomer was considered in a different protonation state as suggested in (1), and in each run the system was simulated for at least 50 ns, using GROMACS 4.0.3 and the GROMOS96 53a6 forcefield. We report three possible proton conduction pathways through the trans-membrane domain. The pathways were identified based on the dynamics of protein-internal water molecules and monitoring their frequency of hydrogen bond formation to adjacent residues. That way we also determined residues likely involved in AcrB's hydrogen bonded network. Each residue was further characterized by its specific hydrogen bond frequency to protein-internal water. Additionally a new method was applied to analyze AcrB-internal cavities and transport tunnels in each monomer's porter domain.

GROMOS96 43a1 performance in predicting oligosaccharide conformational ensembles within glycoproteins

In previous work [Pol-Fachin, L.; Fernandes, C. L.; Verli, H.; Carbohydr. Res. 2009, 344, 491–500], we had demonstrated that GROMOS96 43a1 force field and Löwdin HF/6-31G∗∗-derived atomic charges, adequately represent a glycoprotein’s conformational ensemble in aqueous solutions, taking as the starting geometries NMR-determined structures. Based on such data, the present work intends to evaluate the use of the main solution conformations of isolated disaccharides, to build the carbohydrate moiety of glycoproteins, for which no previous experimental information is available. The observed results suggested that the entire glycoprotein scaffold appears unable to promote major modifications in the conformational behavior of glycosidic linkages. Additionally, when compared to energy contour plots, the results support the use of solution ensembles, to refine vacuum conformations of carbohydrate databases in the assembling of glycoproteins 3D structures. Finally, such approach is applied to build a full glycosylated model for COX-1 and COX-2 enzymes.

Modelling of β- d-glucopyranose ring distortion in different force fields: a metadynamics study

Modelling of carbohydrate conformations is a challenging task for force field developers. Three carbohydrate force fields, namely GLYCAM06, GROMOS 45a4 and OPLS were evaluated. Free energies of different ring conformations of β- d-glucopyranose were calculated using metadynamics in vacuum as well as in explicitly modelled water. All three force fields model the 4 C 1 conformation as the most stable by at least 6 kJ/mol, as compared to other conformations. Interconversion from the 4 C 1 to any other conformation is associated with a barrier of no lower than 26 kJ/mol. The free energy surface calculated in the GLYCAM06 force field is in remarkably good agreement with the recent Car-Parrinello metadynamics study. The effect of a water environment is relatively low and analogous in all tested force fields. Namely, the presence of water stabilizes the upper-left ( 3,O B) versus bottom-right ( B 3,O) area of Stoddard’s plot, relative to the situation in vacuum. Comparison of free and potential surfaces is also provided for vacuum calculations.

Psalmotoxin-1 Docking to Human Acid-sensing Ion Channel-1

Acid-sensing ion channel-1 (ASIC-1) is a proton-gated ion channel implicated in nociception and neuronal death during ischemia. Recently the first crystal structure of a chicken ASIC was obtained. Expanding upon this work, homology models of the human ASICs were constructed and evaluated. Energy-minimized structures were tested for validity by in silico docking of the models to psalmotoxin-1, which potently inhibits ASIC-1 and not other members of the family. The data are consistent with prior radioligand binding and functional assays while also explaining the selectivity of PcTX-1 for homomeric hASIC-1a. Binding energy calculations suggest that the toxin and channel create a complex that is more stable than the channel alone. The binding is dominated by the coulombic contributions, which account for why the toxin-channel interaction is not observed at low pH. The computational data were experimentally verified with single channel and whole-cell electrophysiological studies. These validated models should allow for the rational design of specific and potent peptidomimetic compounds that may be useful for the treatment of pain or ischemic stroke.

Lipid Models for United-Atom Molecular Dynamics Simulations of Proteins.

United-atom force fields for molecular dynamics (MD) simulations provide a higher computational efficiency, especially in lipid membrane simulations, with little sacrifice in accuracy, when compared to all-atom force fields. Excellent united-atom lipid models are available, but in combination with depreciated protein force fields. In this work, a united-atom model of the lipid 1,2-dipalmitoyl-sn-glycero-3-phosphocholine has been built with standard parameters of the force field GROMOS96 53a6 that reproduces the experimental area per lipid of a lipid bilayer within 3% accuracy to a value of 0.623 ± 0.011 nm(2) without the assumption of a constant surface area or the inclusion of surface pressure. In addition, the lateral self-diffusion constant and deuterium order parameters of the acyl chains are in agreement with experimental data. Furthermore, models for 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol (POPG) result in areas per lipid of 0.625 nm(2) (DMPC), 0.693 nm(2) (POPC), and 0.700 nm(2) (POPG) from 40 ns MD simulations. Experimental lateral self-diffusion coefficients are reproduced satisfactorily by the simulation. The lipid models can form the basis for molecular dynamics simulations of membrane proteins with current and future versions of united-atom protein force fields.

GROMOS96 43a1 performance on the characterization of glycoprotein conformational ensembles through molecular dynamics simulations

Considering the small number of papers assessing the conformational profile of glycoproteins through molecular dynamics (MD) simulations, the current work reports on a systematic analysis of the performance of the GROMOS96 43a1 force field and Löwdin HF/6-31G ∗∗-derived atomic charges in the conformational description of glycoproteins. The results substantiate the accuracy of the computational representation of glycoprotein conformational ensembles in aqueous solution based on their agreement to available experimental information, supporting further contributions of computational techniques, mainly MD, in future studies on the characterization of glycoprotein structure and function.

Structural stability of electrosprayed proteins: temperature and hydration effects

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

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.

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, 3 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. 3 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.