Testing of the GROMOS Force-Field Parameter Set 54A8: Structural Properties of Electrolyte Solutions, Lipid Bilayers, and Proteins.

The GROMOS 54A8 force field [Reif et al. J. Chem. Theory Comput.2012, 8, 3705-3723] is the first of its kind to contain nonbonded parameters for charged amino acid side chains that are derived in a rigorously thermodynamic fashion, namely a calibration against single-ion hydration free energies. Considering charged moieties in solution, the most decisive signature of the GROMOS 54A8 force field in comparison to its predecessor 54A7 can probably be found in the thermodynamic equilibrium between salt-bridged ion pair formation and hydration. Possible shifts in this equilibrium might crucially affect the properties of electrolyte solutions or/and the stability of (bio)molecules. It is therefore important to investigate the consequences of the altered description of charged oligoatomic species in the GROMOS 54A8 force field. The present study focuses on examining the ability of the GROMOS 54A8 force field to accurately model the structural properties of electrolyte solutions, lipid bilayers, and proteins. It is found that (i) aqueous electrolytes involving oligoatomic species (sodium acetate, methylammonium chloride, guanidinium chloride) reproduce experimental salt activity derivatives for concentrations up to 1.0 m (1.0-molal) very well, and good agreement between simulated and experimental data is also reached for sodium acetate and methylammonium chloride at 2.0 m concentration, while not even qualitative agreement is found for sodium chloride throughout the whole range of examined concentrations, indicating a failure of the GROMOS 54A7 and 54A8 force-field parameter sets to correctly account for the balance between ion-ion and ion-water binding propensities of sodium and chloride ions; (ii) the GROMOS 54A8 force field reproduces the liquid crystalline-like phase of a hydrated DPPC bilayer at a pressure of 1 bar and a temperature of 323 K, the area per lipid being in agreement with experimental data, whereas other structural properties (volume per lipid, bilayer thickness) appear underestimated; (iii) the secondary structure of a range of different proteins simulated with the GROMOS 54A8 force field at pH 7 is maintained and compatible with experimental NMR data, while, as also observed for the GROMOS 54A7 force field, α-helices are slightly overstabilized with respect to 3(10)-helices; (iv) with the GROMOS 54A8 force field, the side chains of arginine, lysine, aspartate, and glutamate residues appear slightly more hydrated and present a slight excess of oppositely-charged solution components in their vicinity, whereas salt-bridge formation properties between charged residues at the protein surface, as assessed by probability distributions of interionic distances, are largely equivalent in the GROMOS 54A7 and 54A8 force-field parameter sets.

Introduction

(Bio)molecular systems at nonzero temperature are nowadays routinely described as ensembles of structures statistically representative of the phase space sampled by the system under given thermodynamic boundary conditions.[1−3] If the (bio)molecular system of interest is considered at room temperature and in the condensed phase, quantum effects are usually negligible, and the sampling of configurations is thus required to satisfy the Boltzmann distribution. The generation of corresponding configurations can be performed numerically, e.g., based on a molecular dynamics (MD) integration scheme. In classical atomistic MD simulation, the atoms forming the molecular system are represented as mass-possessing point particles associated with a unique interaction site, and an empirical classical potential-energy function (force field) is used to derive the forces acting on the atoms.[4−6] The results obtained from MD simulations depend on the extent of sampling of the phase space accessible to the molecular system and the accuracy of the underlying force field. A major effort in the reign of molecular simulation is thus directed toward the parametrization and validation of force fields against available experimental data.

A coarse distinction of available force fields can be made based on the purpose of the force field, which determines the underlying parametrization strategy. Spectroscopic force fields[7−13] are designed to investigate gas-phase molecular properties, and their parametrization aims at the accurate representation of molecular geometric, energetic, and vibrational properties, while thermodynamic force fields[14−22] are designed to investigate condensed-phase properties, the parametrization accordingly aiming at the accurate representation of bulk thermodynamic properties such as, e.g., densities, vaporization enthalpies, or solvation free energies. The present study is concerned with the GROMOS force field for (bio)molecular simulation. For a thorough discussion of its history and philosophy, the reader is referred to ref (20). The most recent version of the GROMOS force field is the 54A8 parameter set,[22] which is based on (i) the simple-point-charges (SPC) water model;[23] (ii) neutral amino acid side chain parameters from the 53A6 force field,[24] obtained from optimization against pure-liquid properties and solvation free energies of small molecules; (iii) charged amino side chain parameters obtained from optimization against single-ion hydration free energies based on a standard absolute intrinsic proton hydration free energy of[25] ΔGhyd⊖[Hg+] = −1100 kJ·mol–1; (iv) carbohydrate parameters from the 45A4 force field, along with atomic partial charges from ref (26) obtained from restrained electrostatic potential (RESP) fitting compatible with the 53A6 force field; (v) nucleic-acid parameters from the 45A4 force field, along with dihedral angles obtained from quantum-mechanical calculations and atomic partial charges compatible with the 53A6 force field, by ref (27); (vi) lipid parameters from ref (28), along with atomic partial charges from ref (29); (vii) sodium and chloride ion parameters from ref (30) obtained from optimization against single-ion hydration free energies based on a standard absolute intrinsic proton hydration free energy of[25] ΔGhyd⊖[Hg+] = −1100 kJ·mol–1. Besides, alternative parameter sets compatible with the GROMOS force field are available, as e.g. for saturated and unsaturated lipid molecules with phosphocholine or phosphoglycerol head groups,[31] and for hexapyranose-based carbohydrate molecules.[32−34]

Three putative shortcomings of the GROMOS 54A8 force-field parameter set require careful investigation:

1. A crucial feature of the GROMOS 54A8 force-field parametrization strategy is the usage of methodology-independent ion hydration free energies, i.e. the molecular ions analogous to charged amino acid side chain functionalities (guanidinium for arginine, imidazolium for histidine, methylammonium for lysine, and acetate for aspartate and glutamate) were calibrated for the ideal situation of Coulombic electrostatic interactions in a macroscopic nonperiodic system. However, currently, limitations in computing power prohibit the conduction of simulations according to this ideal situation. Typically, systems of finite size and approximate electrostatic interaction functions are resorted to, and the resulting discrepancy between the representation of the charged functional groups during the parametrization procedure (post-simulation correction of raw hydration free energies obtained from configurations sampled with spurious forces for approximate electrostatics, finite-size, and potential-summation artifacts) and during configurational sampling in actual (bio)molecular simulations (no correction for approximate electrostatics, finite-size, and potential-summation artifacts at the level of forces) might lead to artifacts in the configurational sampling. For instance, since the correction terms coining methodology-independent ion hydration free energies are predominantly negative, ions might in effect be underhydrated during simulations, this effect being more pronounced for cations than for anions and presumptively also more pronounced for monatomic ions (concentrated charge density) compared to oligoatomic ones (more diffuse charge density),[22] and a possible consequence of the decreased effective hydrophilicity of charged functional groups might manifest in an overestimation of ion–ion binding in solution, or an underestimation of the interaction of protein side chains with the solvent.

2. Not all charged functional groups were reparameterized in the development of the GROMOS 54A8 force-field parameter set. For instance, the phosphate ion has unaltered Lennard-Jones parameters and atomic partial charges which might not yield a hydration free energy compatible with the above choice of proton hydration free energy and consequently imply imbalances in the relative interaction energies of this ion with other charged functional groups and with water as a solvent.

3. The GROMOS force field relies on a geometric-mean combination rule[35] to define heteroatomic Lennard-Jones interaction parameters based on homoatomic ones. Although a physical rationalization of this rule based on atomic polarizability products in the Slater–Kirkwood formula for dispersion coefficients[36] may be conceived of in the case of interactions between isolated monatomic species, it remains truly ad hoc for repulsion coefficients and in the case of interactions between atoms in molecules, and one is hence bound to question its tenability upon translation of Lennard-Jones interaction parameters valid in a system governed by few heteroatomic interactions (e.g., a single acetate ion in aqueous solution) to an environment characterized by multiple interactions with molecular species and functional groups other than water [e.g., an acetate ion in a 1.0 m (1.0-molal) aqueous sodium acetate solution, or a carboxylate group pertaining to a glutamate residue in a protein].

In the context of (bio)molecular simulation, it has therefore to be tested whether, e.g., biomacromolecular structure and dynamics, the structural, dynamic, and thermodynamic properties of electrolyte solutions, and interactions between the two latter types of systems as, e.g., encountered in the province of Hofmeister effects[37−40] on biomacromolecules, are accounted for in a realistic way. Considering aqueous electrolyte solutions involving ionic species parametrized against methodology-independent hydration free energies, a realistic description of the thermodynamic equilibrium between ion–ion and ion–water pairing propensities in the GROMOS 54A8 force field is not a priori guaranteed because (i) interionic distances and water configurations sampled in systems of finite sizes and with approximate electrostatic interaction functions might lead to artificial interionic distance distributions and inadequate ion hydration[41]and (ii) the geometric-mean combination rule applied on the basis of optimal ion–water Lennard-Jones interactions might prove inadequate for the representation of ion–ion Lennard-Jones interactions.[42]

Considering arrangements of lipid molecules in solution, the increased van der Waals radius of the methyl group attached to positively charged species is expected to be reflected in structural parameters such as, e.g., area per lipid of choline-headgroup containing lipids. In the GROMOS 54A7 force field,[21] the repulsion between choline methyl groups and phosphate oxygen atoms was increased to reproduce experimental properties for the liquid crystalline-like phase.[28] In the GROMOS 54A8 force field, this change was reverted, because the increased Lennard-Jones repulsion parameter of the choline methyl group in comparison to the GROMOS 54A7 force field achieves approximately the same alteration in lipid headgroup repulsion (section III.2).[22]

Considering native proteins in solution, the most decisive signature of the GROMOS 54A8 force field in comparison to its predecessor 54A7 can probably be found in the thermodynamic equilibrium between salt-bridge formation and hydration of charged functional groups at the protein surface. Since salt bridges influence protein stability only marginally,[43−46] it can be expected that the GROMOS 54A8 force field yields stable proteins akin to the GROMOS 54A7 force field. Similar considerations apply to the situation of proteins in aqueous electrolyte solutions. On the basis of the law of “matching water affinities”,[40,47,48] it may be presumed that weakly hydrated anionic species (e.g., chloride ions) bind directly to positively charged functional groups in a protein (which are likewise weakly hydrated) and that strongly hydrated cationic species (e.g., sodium ions) bind directly to protein carboxylate groups (which are likewise strongly hydrated). However, failure to properly account for the underlying free energies of hydration may confound the relative ion pairing propensities, which could entail an erroneous classification of ions concerning their effect on protein stability (“Hofmeister series,”[38] e.g. salting-in vs salting-out phenomena).

The goal of this article is to investigate some of the above issues in simulations of electrolyte solutions, lipid bilayers, and proteins. In particular, it is studied whether (i) thermodynamic properties (experimental activity derivatives) of electrolyte solutions are correctly reproduced; (ii) usage of the GROMOS 54A8 force-field representation of the tetramethylammonium ion[22] in the positively charged choline headgroup of phosphatidylcholine lipids still achieves a realistic description of lipid bilayer structural properties; (iii) the GROMOS 54A8 force field yields stable protein secondary structures compatible with experimental NMR data (NOE values and 3J-coupling constants; where available); (iv) protein surface salt-bridge formation and counterion binding properties differ between the GROMOS 54A7 and 54A8 force fields. The results of these investigations are expected to be sensitive to the employed simulation methodology. In order for them to be of highest relevance to the GROMOS user community, the most widely used simulation protocol is employed,[49] which implies, notably, a cutoff-truncated electrostatic interaction function with reaction-field correction. However, in the case of simulations of electrolyte solutions, an exception to the standard protocol is made, and a lattice-sum electrostatic interaction function is used. This was necessary to allow a valid comparison with experimental data for those highly electrostatic artifact-prone systems.

Computational Details

Setup of Simulated Systems

The simulated systems and performed simulations are summarized in Table 1.

Table 1

Abbreviations for the Simulated Systems and Performed Simulations Used Throughout the Texta

abbreviation	system and simulation specification	force field	simulation time [ns]	⟨L⟩ [nm]
electrolyte solutions
NA-CL	0.5, 1.0, or 2.0 m aqueous sodium chloride solution	54A8	6	5.98,b 5.95,c 5.93d
NA-ACET	0.5, 1.0, or 2.0 m aqueous sodium acetate solution	54A8	6	6.00,b 5.98,c 5.98d
H3C1-CL	0.5, 1.0, or 2.0 m aqueous methylammonium chloride solution	54A8	6	6.02,b 6.03,c 6.07d
GUAN-CL	0.5, 1.0, or 2.0 m aqueous guanidinium chloride solution	54A8	6	6.03,b 6.04,c 6.10d
lipid bilayers
DPPC1	hydrated DPPC bilayer	54A8	40	6.41,e 8.29f
DPPC2	hydrated DPPC bilayer	g	40	6.62,e 7.79f
DPPC3	hydrated DPPC bilayer	h	40	6.73,e 7.53f
proteins
HEWL	hen egg white lysozyme in water	54A8	20	7.72
FOX	RNA-binding domain of the human FOX-1 protein in complex with a UGCAUGU RNA strand in water	54A8	20	7.18
CM	chorismate mutase of Mycobacterium tuberculosis in water	54A8	20	8.86
GCN	GCN4 trigger peptide in water	54A8	100	4.48
PROTG	B1 immunoglobulin-binding domain of streptococcal protein G in water	54A8	20	6.19
COLDS	major cold chock protein CspA of Escherichia coli in water	54A8	20	6.38
SAC054A7	hyperthermophilic protein Sac7d of Sulfolobus acidocaldarius (SAC) in water	54A7	20	6.60
SAC054A8	SAC in water	54A8	20	6.60
SACCl–54A7	SAC in water, with six chloride counterions	54A7	20	6.60
SACCl–54A8	SAC in water, with six chloride counterions	54A8	20	6.60
SACNaCl54A7	SAC in water, with six sodium and 12 chloride counterions	54A7	20	6.60
SACNaCl54A8	SAC in water, with six sodium and 12 chloride counterions	54A8	20	6.60

The indicated simulation time refers to the length of the production run, and the reported box-edge lengths ⟨L⟩ are arithmetic averages over the production run. If not indicated otherwise, a neutralizing amount of sodium or chloride counterions was used in the protein simulations. More details are provided in sections II.1.1, II.1.2, and II.1.3 for simulations involving electrolyte solutions, lipid bilayers, and proteins, respectively.

0.5 m.

1.0 m.

2.0 m.

Box-edge length along the bilayer plane.

Box-edge length along the bilayer normal.

Modified GROMOS 54A8 force-field parameter set, where the C12 parameter of atom type 2 is employed for interactions with atom type 54.

Modified GROMOS 54A8 force-field parameter set, where the C12 parameter of atom type 2 is employed for interactions with atom type 54.

Electrolyte Solutions

Aqueous 0.5, 1.0, and 2.0 m solutions of sodium chloride (NA-CL), sodium acetate (NA-ACET), methylammonium chloride (H3C1-CL), and guanidinium chloride (GUAN-CL) were created by randomly placing Ns ion pairs and Nw water molecules (Table 2), the latter represented by the SPC water model,[23] in cubic computational boxes whose sizes were as a first approximation chosen such that the overall molecule number density satisfied the equilibrated number density of the SPC water model (32.5 nm–3).[50] An equilibration run consisting of 4 ns with cutoff-truncated electrostatic interactions including a Barker–Watts reaction-field correction and 1 ns with LS electrostatic interactions was performed in the NPT ensemble at a pressure of 1 bar and a temperature of 298.15 K. Further details concerning the handling of electrostatic interactions are provided in section II.2.

Table 2

Computed and Experimental Activity Derivatives fss (eq 2) for the Salts NA-CL, NA-ACET, H3C1-CL, and GUAN-CL Corresponding to Aqueous Electrolyte Solutions of Salt Molality bsa

	bs [m]	xs	Ns	Nw	Gss [nm3]	Gsw [nm3]	Gww [nm3]	⟨L3⟩ [nm3]	fss (sim.)	fss (exp.)
NA-CL	0.5	0.0177	62	6881	0.64	–0.013	–0.029	214.0	–0.27	–0.055
NA-ACET					–0.11	–0.025	–0.029	215.7	0.05	0.017
H3C1-CL					0.011	–0.044	–0.028	218.5	–0.04	–0.106
GUAN-CL					0.21	–0.060	–0.028	219.5	–0.14	–0.194
NA-CL	1.0	0.0348	121	6715	0.66	–0.026	–0.029	210.3	–0.43	0.011
NA-ACET					–0.11	–0.023	–0.029	213.5	0.11	0.135
H3C1-CL					–0.038	–0.043	–0.028	218.9	–0.02	–0.063
GUAN-CL					0.12	–0.063	–0.027	220.8	–0.19	–0.246
NA-CL	2.0	0.0672	235	6521	1.69	–0.092	–0.026	208.1	–0.80	0.169
NA-ACET					–0.12	–0.020	–0.029	213.7	0.29	0.354
H3C1-CL					–0.090	–0.038	–0.028	224.0	0.09	0.046
GUAN-CL					–0.008	–0.062	–0.024	227.5	–0.15	–0.302

The reported data include the number of cations Ns, the number of anions Ns, the number of water molecules Nw, the corresponding salt mole fraction xs (eq A.3), the average volume of the computational box ⟨L3⟩, as well as the Kirkwood–Buff integrals G for salt–salt (ss), salt–water (sw), and water–water (ww) interactions. The experimental activity derivatives are obtained from the fitted functions given by eq A.2 using the coefficients a, b, c, and d reported in the Supporting Information, Table SI. Simulated activity derivatives are calculated according to eq 5, using Kirkwood–Buff integrals G given by the y-axis intercept of a linear fit to G(Ls–1), as displayed in Figure 1.

Figure 1

Kirkwood–Buff integrals G(Ls–1) for salt–salt (ss), salt–water (sw), and water–water (ww) interactions calculated from MD simulations of the systems listed in Table 2 for salt molalities bs = 0.5 m (black), 1.0 m (red), and 2.0 m (green). Circles represent G(Ls–1) evaluated via eq 4 using sub-boxes of edge length Ls randomly positioned in computational boxes of average volume ⟨L3⟩ reported in Table 2. Solid lines represent a linear fit to the linear regime of G(Ls–1). The choice of linear regime appears ambiguous for salt–salt and salt–water interactions in 2 m NA-CL, which is why for these systems the fit was performed over the range Ls–1 > Ls,max–1, where Ls,max–1 denotes the onset of divergent behavior in G(Ls–1), i.e., approximately 0.39 nm–1 (section III.1).

Lipid Bilayers

A configuration of a hydrated dipalmitoylphosphatidylcholine (DPPC) bilayer obtained after 190 ns of MD simulation with the GROMOS 54A7 force field at a pressure of 1 bar and a temperature of 323 K by ref (28) was downloaded from the automated topology builder repository.[51] It contains 128 DPPC molecules (8 × 8 molecules per leaflet) in the liquid crystalline-like phase and 5841 SPC water molecules in a box of dimensions 6.3398 × 6.3398 × 8.3985 nm3, the bilayer plane being parallel to the xy plane of the computational box. The area per lipid AL (eq 7) of this configuration thus evaluates to 0.628 nm2, close to the experimental values reported by ref (52) (0.629 nm2) and ref (53) (0.630 nm2) and somewhat lower than the values reported by ref (54) (0.64 nm2) and ref (55) (0.642–0.643 nm2). This configuration was used as the starting configuration in three MD simulations denoted DPPC1, DPPC2, and DPPC3, which differ in the C12 parameter of the phosphate oxygen atom type (atom type 2 in the GROMOS 54A7 and 54A8 force fields) to be used for interactions with the methyl group in the choline headgroup (atom type 54 in the GROMOS 54A7 and 54A8 force fields), namely [C12]1/2 = 0.8611 × 10–3, [C12]1/2 = 1.841 × 10–3, or [C12]1/2 = 3.068 × 10–3 kJ1/2·mol–1/2·nm6, respectively. The first choice was adopted in the construction of the GROMOS 54A8 force field, i.e., simulation DPPC1 corresponds to this set of force-field parameters (section I), whereas simulations DPPC2 and DPPC3 present increased Lennard-Jones repulsion between the choline head groups and the phosphate oxygen atoms. The third choice (albeit in combination with a C12 parameter of atom type 54 that is lower than in the GROMOS 54A8 force field) was adopted in the construction of the GROMOS 54A7 force field.[21,28] Note, in addition, that unlike previous GROMOS force-field parameter sets, the GROMOS 54A7 and 54A8 force-field representations of the DPPC molecule adopt the atom partial charge distribution of Chiu et al.[29] These charges are derived from ab initio quantum-mechanical calculations and turn out to be significantly enhanced in magnitude in comparison to their GROMOS pre-54A7 analogs. Although adherence to quantum-mechanically derived charge distributions is actually at variance with the philosophy of the GROMOS force field,[3,35] it turned out necessary in this specific case to adequately model the biologically important liquid crystalline-like phase of DPPC bilayers.[28]

The structure downloaded from the automated topology builder repository[51] corresponds to a temperature of 323 K, but as no velocities were available it was equilibrated in six successive MD runs each of 20 ps length in the NVT ensemble, for each of the DPPC1, DPPC2, and DPPC3 simulations. The three simulations were initiated from a random set of atom velocities satisfying a Maxwell–Boltzmann distribution at 60 K. In the first simulation period, harmonic position restraints on the lipid atoms were applied with a force constant of 2.5 × 104 kJ·mol–1·nm–2. In the subsequent simulations, the temperature was increased to 120, 180, 240, 300, and 323 K while the force constant was reduced to 2.5 × 103, 2.5 × 102, 25, 2.5, and 0 kJ·mol–1·nm–2, respectively.

Proteins

Hen egg white lysozyme (HEWL), the RNA-binding domain of the human FOX-1 protein in complex with a UGCAUGU RNA strand (FOX), chorismate mutase of Mycobacterium tuberculosis (CM), the GCN4 trigger peptide (GCN), the B1 immunoglobulin-binding domain of streptococcal protein G (PROTG), the major cold chock protein CspA of Escherichia coli (COLDS), and the hyperthermophilic protein Sac7d of Sulfolobus acidocaldarius (SAC) were simulated with the GROMOS 54A8 force field to investigate secondary structure stability and, where available, compatibility of simulation results with experimental NMR data. SAC was in addition simulated with the GROMOS 54A7 force field because it contains a multitude of charged surface residues. For this protein, a comparison between both force fields in terms of surface salt-bridge formation and counterion binding properties thus appears interesting.

Starting configurations for HEWL, FOX, CM and GCN were provided by the authors of ref (21) who used these systems for testing the GROMOS 54A7 force field, i.e., these starting configurations are rigorously identical to those used in ref (21) and correspond to entries 1AKI (HEWL; X-ray diffraction structure at 1.5 Å resolution),[56]2ERR (FOX; entry number 1 from a solution-NMR structure bundle),[57]2FP2 (CM; X-ray diffraction structure at 1.6 Å resolution),[58] and 2OVN (GCN; entry number 1 from a solution-NMR structure bundle)[59] in the Protein Data Bank (PDB). Starting configurations for COLDS and PROTG were provided by a collaborator in GROMOS force-field development (Andreas P. Eichenberger, personal communication) and correspond to PDB entries 1MJC (COLDS; X-ray diffraction structure at 2.0 Å resolution)[60] and 1PGB (PROTG; X-ray diffraction structure at 1.9 Å resolution).[61] Finally, a starting configuration for SAC was derived from PDB entry 1SAP (solution-NMR structure).[62] This PDB entry does not provide experimental Nuclear Overhauser Enhancement (NOE) values for SAC. In the following, the details of the system preparation protocol for 1SAP are provided and are, for the sake of completeness, also briefly repeated for the other proteins.

The X-ray diffraction structures (in the case of CM, only the first subunit of the reported dimer) and the first structures of the NMR bundles (FOX, GCN) were complemented with missing hydrogen coordinates if necessary, relaxed in a vacuum using the steepest-descent energy minimization functionality of the GROMOS11 program,[49] and solvated in cubic computational boxes composed of periodically replicated cubic water boxes containing 216 SPC water molecules at the equilibrated density of the SPC water model (972 kg·m–3),[50] allowing for a minimum protein–solvent distance of 0.23 nm. The resulting box-edge lengths are reported in Table 1. To remove steric strain involving the solvent molecules, another steepest-descent energy minimization was applied. During this procedure, the solute atom positions were restrained to their initial positions using a force constant 2.5 × 104 kJ·mol–1·nm–2. In simulations HEWL, CM, PROTG, and COLDS, a neutralizing amount of sodium or chloride counterions was added by replacing randomly chosen water molecules within the bulk solvent. In simulations FOX and GCN, no counterions were added to maintain consistency with previous studies involving these proteins.[21,63] All simulations were performed with amino acid side chain and terminal-group protonation states appropriate for pH 7, i.e., all arginine, lysine, aspartate, and glutamate residues, as well as protein termini were charged. One exception is formed by residue K5 of the FOX protein, which was kept deprotonated to maintain consistency with previous studies involving this protein[21,63] where the deprotonation of this residue was erroneously introduced. The resulting number of charged side chains is reported in Table 4, and it follows that simulations HEWL, FOX, CM, GCN, PROTG, and COLDS involve net charges of 0, −4, 0, +1, 0, and 0e, respectively, where the charge of −6e of the RNA molecule binding to the FOX protein was taken into account. For SAC, the system preparation followed the same protocol, except for the following: (i) systems were prepared for both GROMOS 54A7 and 54A8 force fields (corresponding simulations are denoted with superscripts “54A7” and “54A8”; Table 1); (ii) the protein was not minimized in vacuo to avoid possible structural disruptions which might arise as a consequence of spurious salt-bridge formation between the many charged surface residues. Instead, a steepest-descent energy minimization was done for the entire system after solvation; (iii) addition of counterions was not only done with the aim of neutralizing the protein, which carries a net charge of +6e. Instead, for each force-field parameter set, three approaches were investigated, namely omission of counterions (denoted SAC0), addition of a neutralizing amount of chloride ions (denoted SACCl), or addition of a neutralizing sodium chloride solution consisting of 12 chloride and six sodium ions (denoted SACNaCl).

Table 4

Parameters Characterizing Stability of Secondary Structure Evaluated from the Indicated Proteinsa

simulation	rmsdbb (max.) [nm]	Nbb–hb (% of ini.)	fα-helix [%]	f310-helix [%]	fπ-helix [%]	fβ-sheet [%]	NR	NK	ND	NE
HEWL	0.20 (0.31)	59.1 (89.5)	34.4	3.7	1.7	6.1	11	6	7	2
ini.		66	29.5	14.0	0.8	6.2
FOX	0.29b (0.39)b	34.9 (79.3)	17.5	1.5	0.0	21.8	9	5 (6)c	5	7
ini.		44	21.6	0.0	0.0	23.9
CM	0.33 (0.44)	102.3 (87.4)	74.1	2.3	0.0	0.0	14	4	14	9
ini.		117	72.4	6.7	0.0	0.0
GCN	0.34 (0.55)	8.0 (80.0)	47.3	0.3	16.6	0.0	1	2	0	2
ini.		10	62.5	12.5	0.0	0.0
PROTG	0.13 (0.32)	29.5 (95.2)	25.9	0.1	0.1	39.1	0	6	5	5
ini.		31	25.0	0.0	0.0	42.9
COLDS	0.17 (0.29)	30.1 (88.5)	0.0	3.3	0.0	42.1	0	7	6	2
ini.		34	0.0	4.3	0.0	40.6
SAC054A7	0.27 (0.37)	28.9 (80.3)	16.0	3.7	0.1	35.7	4	14	5	7
SAC054A8	0.27 (0.50)	28.9 (80.3)	16.7	3.2	0.0	37.0
SACCl–54A7	0.22 (0.40)	27.7 (76.9)	16.7	0.6	0.0	34.7
SACCl–54A8	0.28 (0.45)	28.5 (79.2)	18.9	2.5	0.0	31.6
SACNaCl54A7	0.23 (0.34)	29.4 (81.7)	17.9	2.4	0.1	36.0
SACNaCl54A8	0.36 (0.52)	25.8 (71.7)	13.2	2.8	0.0	34.4
ini.		36d	21.2d	9.1d	0.0d	45.5d

The reported data include the average backbone heavy atom–positional rmsd from the structure obtained after equilibration, rmsdbb, evaluated after fitting the backbone Cα atoms to this structure, as well as the maximum instantaneous rmsd value in parentheses, the average number of backbone hydrogen bonds Nbb–hb, as well as the corresponding fraction of the initial number of backbone hydrogen bonds in parentheses, and average fractions fα-helix, f3, fπ-helix, and fβ-sheet of protein residues occurring in secondary structure elements α-helix, 310-helix, π-helix, or β-sheet, respectively. Table entries “ini.” denote data corresponding to the initial (minimized) structures. In addition, the number of charged arginine (NR), lysine (NK), aspartate (ND), and glutamate (NE) side chains is reported.

Considering residues 9–82 only (the rmsd of the backbone atoms of all residues, i.e. residues 1–88, is significantly higher due to large fluctuations in the tail regions and its average and maximum values evaluate to 0.54 and 0.68 nm, respectively).

The FOX protein contains six lysine residues, of which only five were protonated to maintain consistency with previous studies involving this protein,[21,63] where the deprotonation of this residue was erroneously introduced (section II.1.3).

This number turns out to be the same regardless of whether the minimization of the system was performed with the GROMOS 54A7 or 54A8 force-field parameter set.

All simulations were initiated from a random set of atom velocities satisfying a Maxwell–Boltzmann distribution at 60 K, and the systems were equilibrated in five successive MD runs each of 20 ps length in the NVT ensemble. In the first simulation period, harmonic position restraints on the protein atoms were applied with a force constant of 2.5 × 104 kJ·mol–1·nm–2. In the subsequent simulations, the temperature was increased to 120, 180, 240, and 300 K, while the force constant was reduced to 2.5 × 103, 2.5 × 102, 25, and 0 kJ·mol–1·nm–2, respectively. For CM, a temperature of 310 K was used in the fifth run to maintain consistency with a previous study[21] involving this protein.

Molecular Dynamics Simulations

Twelve simulations of electrolyte solutions, three simulations of hydrated lipid bilayers, and 12 simulations of hydrated proteins were performed using the GROMOS11 package of programs[49] (Table 1). Production runs were started from equilibrated system configurations obtained as described in section II.1. They lasted 6 ns for the electrolyte solutions, 40 ns for the lipid bilayers, 20 ns for all protein systems except GCN, and 100 ns for GCN. See Supporting Information section S.I for details on the simulation settings.[35,64−78]

Analysis of Simulation Results

Calculation of Electrolyte Activity Derivatives

The rational (mole-fraction scale) activity coefficient fs of a salt relates the salt activity as to the salt mole fraction xs, i.e.[79]The activity derivative fss is given byand can be evaluated from the structural characteristics of the electrolyte solution by means of a Kirkwood–Buff analysis, because Kirkwood–Buff theory[80−82] establishes the necessary connection between thermodynamic and structural solution properties. The latter are characterized through Kirkwood–Buff integrals[80]G, which may be calculated on the basis of radial distribution functions (rdfs) gμ(r) of particles I and J obtained from MD simulations in the μVT ensemble but are generally approximated using rdfs g(r) corresponding to the NPT ensemble, i.e.or on the basis of particle number fluctuations in the μVT ensemble, i.e.where N and N denote the number of particles I and J, respectively, angular brackets denote an ensemble average, and δ denotes the Dirac delta function (δ = 1 for I = J and zero otherwise).

For the approach based on eq 3, the integration is in practice performed from 0 to finite distances r smaller than half the edge length of the computational box, and the Kirkwood–Buff integrals G(r) = 4π ∫0 dr′[g(r′) – 1]r′2 are averaged over a range of values r for which the rdfs are reasonably close to the bulk value g(∞) = 1. If I = J and the rdfs are obtained from a simulation at constant particle number, a finite-size correction (typically via scaling of the rdf) has to be applied to account for the reduced like-particle density around the central particle.[48]

The approach based on eq 4 was suggested[83,84] to be performed in a sub-box of edge length Ls cut out from a larger box for which a simulation trajectory pertaining to a constant-particle-number ensemble is available. In this case, the volume V in eq 4 is to be replaced by the volume of the sub-box, Ls3, and the particle numbers N and N in eq 4 are counted in the sub-box. Thus, N and N are not constant, i.e., effectively correspond to an open ensemble. Finite-size effects may be accounted for by evaluating eq 4 for multiple values of Ls and extrapolating the linear regime of G(Ls–1) to Ls–1 = 0.[83,84] The authors of the present study experienced difficulties with the convergence of G(r) calculated via eq 3 for salt–water and salt–salt interactions. Therefore, G(Ls–1) was calculated via eq 4 for randomly positioned cubic sub-boxes of edge length Ls ranging from 0.5 to 1.0 nm in steps of 0.05 nm and from 1.0 to 5.25 nm in steps of 0.25 nm, and a finite-size corrected value G was obtained as the y-axis intercept of a linear fit to G(Ls–1). Although small values of Ls allow analysis of multiple nonoverlapping sub-boxes per trajectory frame, this was not done in the present study because plain analysis using one sub-box per trajectory frame already provided sufficient statistics. The counting of oligoatomic ions and water molecules in the sub-box relied on the carboxylate carbon atoms of acetate ions, the nitrogen atoms of methylammonium ions, the carbon atoms of guanidinium ions, and the oxygen atoms of water molecules. For most of the salts and concentrations examined, the resulting Kirkwood–Buff integrals turned out to be similar to the values obtained with the rdf-based approach (eq 3) if the above-mentioned averaging of G(r) was done over the interval 1.0 ≤ r ≤ 1.5 nm (data not shown).

The activity derivative fss is related to the salt–salt (ss), salt–water (sw), and water–water (ww) Kirkwood–Buff integrals as[81]where ηw is the number of density of water molecules andNote that the above outline of the application of Kirkwood–Buff theory to an electrolyte solution treats cations and anions as indistinguishable particles.[85,86] This approach allows a meaningful comparison to experimental activity derivative data, the corresponding activity coefficients being, due to the electroneutrality of macroscopic matter at equilibrium, mean rather than single-ion values, i.e. fs2 = ff for a 1:1 electrolyte composed of ions Iaq and Jaq–.

Experimental data for fss can be obtained by differentiation of fitted functions for experimental fs values (Appendix, eq A.7).

Biomolecular Structure

The structural properties of the simulated DPPC lipid bilayers were analyzed in terms of area and volume per lipid and electron density profiles. The average area per lipid AL and volume per lipid VL were calculated asandrespectively, where NL = 128 is the number of DPPC molecules, LL is the area of one lipid bilayer leaflet in the computational box, calculated as the product of the box-edge lengths along the x and y directions, and VBL = LLL – NwVH2O is the volume of the computational box occupied by the lipid bilayer, calculated as the box volume LLL reduced by an estimate NwVH2O for the volume of the box occupied by water. The latter is approximated as the product of the number of water molecules Nw = 5841 and an estimate for the volume of one SPC water molecule VH2O under the given pressure and temperature conditions (VH2O = 3.15 × 10–2 nm3 at P = 1 bar and T = 323 K).[28] The electron density profile ηe(z) along the bilayer normal direction was calculated aswhere z′ denotes the signed distance of an atom from the bilayer center, ne(P) = 16, ne(O) = 8, ne(N) = 6, and ne(C) = 6 are the numbers of electrons of phosphorus, oxygen, nitrogen, and carbon atoms, respectively, and ⟨NP(z;z′)⟩, ⟨NO(z;z′)⟩, ⟨NN(z;z′)⟩, and ⟨NC(z;z′)⟩ are the average numbers of these atoms sampled at distances |z| – (1/2)Δz ≤ |z′| < |z| + (1/2)Δz from the bilayer center. Electron numbers of 16 and 6 for phosphorus and nitrogen, respectively, account for the net charge of the phosphate and choline groups, respectively. Electrons of hydrogen atoms associated with the united-atom methylene and methyl groups were neglected. z ranges from −L/2 + Δz/2 to L/2 – Δz/2, and the bin width Δz is L/100. The electron density profile of a lipid bilayer in the liquid crystalline-like phase presents two maxima whose positions are mostly determined by the electron-rich phosphorus atoms. The distance between maxima in the electron density of lipid headgroup atoms, ηe(H)(z), is referred to as the head–head bilayer thickness DHH.[28] Here, the headgroup atoms chosen for calculation of the electron density profile involve the choline and the phosphate group, i.e., three united-atom methyl groups, a nitrogen atom, two united-atom methylene groups, one phosphorus atom, and four oxygen atoms, such thatwhere ⟨NO∈H(z;z′)⟩ and ⟨NC∈H(z;z′)⟩ refer to oxygen and carbon atoms in the lipid head groups.

Stability of a protein during the simulation was assessed by (i) monitoring the time series of atom-positional root-mean-square deviation (rmsd) of all backbone heavy atoms (N, Cα, C) from the structure obtained after equilibration after a roto-translational fit of all Cα atoms to this reference structure; (ii) evaluating the time average of the number of backbone hydrogen bonds, based on defining a hydrogen bond via a geometric criterion: a maximum hydrogen atom-acceptor distance of 0.25 nm combined with a minimum donor-hydrogen atom-acceptor angle of 135°; (iii) monitoring the time series of fractions of protein residues occurring in secondary structure elements α-helix, 310-helix, π-helix, or β-sheet, based on secondary structure definitions according to the Kabsch–Sander rules.[87]

For three proteins, the simulated trajectories were also validated against available experimental NMR data, namely, the NOE values of ref (88), ref (57), and refs (59), (89) for HEWL, FOX, and GCN, respectively, and 3J-coupling constants of ref (88) and refs (59), (89) for HEWL and GCN, respectively. To perform a comparison against upper bounds for interproton distances derived from NOE experiments, an inverse-sixth power averaging was applied to the interproton distances rH monitored during the simulation, i.e.where angular brackets indicate an ensemble average. Since the GROMOS force field relies on a united-atom representation of aliphatic CH (n = 1, 2, 3, 4) groups, the positions of hydrogen atoms which are not treated explicitly were deduced based on geometric criteria.[35] Furthermore, in cases where experimental NOE data refer to stereochemically equivalent and nonstereoassigned protons, pseudoatom positions[90] were constructed according to geometric criteria and were used in combination with corresponding standard GROMOS corrections[35] but without any additional multiplicity corrections.[91] Average NOE violations were calculated as the sum of positive violations divided by the number of available NOE values. To perform a comparison against experimental 3J-coupling constants between two protons, these quantities were calculated from the simulated trajectory based on the average dihedral angle φ formed by the respective set of bonded atoms. Here, exclusively coupling constants 3JH(ϕ) between backbone amide hydrogen and Cα hydrogen atoms were considered. However, because the aliphatic Hα, hydrogen atom is merged together with the atom Cα, in the GROMOS force field,[35] the dihedral angle ϕ′ formed by the atoms C–N–Cα,–C was used instead. The 3J-coupling constants were calculated from the Karplus relation[92] aswhere ϕ = ϕ′ – 60° and the parameters of ref (93) were used (A = 6.4 Hz, B = −1.4 Hz, C = 1.9 Hz).

For the simulations of the SAC protein, various pair distribution functions involving charged functional groups were calculated. Radial distribution functions g(r) were calculated for water oxygen atoms, sodium ions, and chloride ions around (i) the amino nitrogen atoms of the alkylguanidinium group of arginine residues, (ii) the nitrogen atoms of the alkylammonium group of lysine residues, and (iii) the carboxylate oxygen atoms of aspartate and glutamate residues, aswhere r denotes the minimum-image distance between groups I and J, η is the number density of group J in the system, ⟨N(r;r)⟩ is the average number of groups J with r – (1/2)Δr ≤ r < r + (1/2)Δr, and the bin width Δr = 0.01 nm. Distance probability distributions Pαβ(r), α,β = +, – denoting the set of charged protein functional groups, were calculated from the distance time series of all pairs of like-charged positive (P++), like-charged negative (P––), and oppositely charged (P+–) protein functional groups. The distance measurements were based on the alkylguanidinium carbon atoms of arginine residues, the alkylammonium nitrogen atoms of lysine residues, and the carboxylate carbon atoms of aspartate and glutamate residues, and the functions Pαβ(r) were normalized such thatwhere the number of distinguishable pairs Npairαβ evaluates to Npair++ = 153, Npair–– = 66, and Npair+– = 216.

All the above analyses were done with the gromos++ package of programs.[94]

Results

Electrolyte Solutions

Average box volumes monitored during the production runs are reported in Table 2. Kirkwood–Buff integrals G(Ls–1) are displayed in Figure 1 for salt–salt, salt–water, and water–water interactions as a function of inverse edge length Ls–1 of the sub-box (eq 4). Table 2 reports the y-axis intercepts of linear fits to G(Ls–1), corresponding to finite-size corrected Kirkwood–Buff integrals. Divergence of G(Ls–1) for Ls–1 → 0 is due to insufficient size of the particle reservoir surrounding the sub-box. Significant deviations from linearity are observed for 0.5 < Ls–1 < 1.0 nm–1 in the case of salt–salt and salt–water Kirkwood–Buff integrals for 2 m NA-CL and 1.5 < Ls–1 < 2.0 nm–1 in the case of salt–water Kirkwood–Buff integrals for 2 m NA-CL, which might be caused by spurious ion clustering in this system. From eq 4, it can be seen that excessive pairing between I and J species (i.e., high correlation between N and N) leads to an increase in G. The increase in the magnitude of the slope of G(Ls–1) for 0.5 < Ls–1 < 1.0 nm–1 in salt–salt and salt–water Kirkwood–Buff integrals for 2 m NA-CL, which leads to an increase of the salt–salt and a decrease of the salt–water Kirkwood–Buff integrals, is thus indicative of increased interionic and decreased ion–water interactions. Similarly, the slope of G(Ls–1) for salt–water interactions in 2 m NA-CL is more negative in the range 1.0 < Ls–1 < 1.5 nm–1 than in the range 1.5 < Ls–1 < 2.0 nm–1, which indicates a corresponding loss of ion–water pairing in sub-boxes of the former size. The choice of linear regime is thus ambiguous for salt–salt and salt–water interactions in 2 m NA-CL. Therefore, for these systems, the fit was performed over the range Ls–1 > Ls,max–1, where Ls,max–1 denotes the onset of divergent behavior in G(Ls–1), i.e., approximately 0.39 nm–1 (Figure 1). Note that the particular choice of fitting range does not crucially affect fss for 2 m NA-CL. For example, a fit in the range 0.5 < Ls–1 < 1.0 nm–1 leads to Gss = 1.97 and Gsw = −0.103 nm3, i.e., to results reflecting enhanced ion–ion pairing at the cost of ion–water pairing in comparison to the results obtained from the fit over the range Ls–1 > Ls,max–1 (Gss = 1.69 and Gsw = −0.092 nm3; Table 2). The fss values corresponding to the former and latter Kirkwood–Buff integrals differ only marginally (−0.82 and −0.80, respectively; Table 2).

Evaluation of the Kirkwood–Buff integrals G for salt–salt, salt–water, and water–water interactions allows calculation of activity derivatives fss (eq 2) via eq 5. These values, along with corresponding experimental data, are reported in Table 2. Activity derivatives fss express the change of the activity coefficient fs with salt concentration. The experimental values of ln fs are displayed in Figure 2a as a function of ln xs. A deviation of ln fs from zero indicates nonideality, which, in a salt solution, may be caused by ion pairing (mutual excess coordination of oppositely-charged salt components) or increased water density in the vicinity of ions (excess coordination of salt components by water molecules). An increase in salt concentration is expected to enhance ion pairing and to reduce hydration effects, i.e., to enhance opposing effects on solution ideality. As illustrated in Figure 2a, an increase in the salt mole fraction has a very limited influence on the (mean rational) activity coefficient up to salt mole fractions xs (eq A.3) of about 0.0177 (ln xs = −4.03; salt molality bs = 0.5 m) for GUAN-CL, of about 0.0348 (ln xs = −3.36; bs = 1.0 m) for NA-ACET and of about 0.0672 (ln xs = −2.70; bs = 2.0 m) for NA-CL and H3C1-CL. The derivatives fss of ln fs with respect to ln xs (eq 2) are displayed in Figure 2b as a function of ln xs for both the simulated and experimental salt solutions, the corresponding numerical values being reported in Table 2. Simulated and experimental data for fss are in very good quantitative agreement for salts NA-ACET, H3C1-CL, and GUAN-CL at concentrations up to 1.0 m (mean and maximum absolute deviations from experiment amounting to 0.05 and 0.07, respectively), and also for NA-ACET and H3C1-CL at 2.0 m concentration (mean and maximum absolute deviations from experiment amounting to 0.05 and 0.06, respectively), whereas a higher deviation from experimental data is found for GUAN-CL at 2.0 m concentration (0.15). The salts involving oligoatomic ionic species thus reproduce experimental data very well, however, not even qualitative agreement is found for NA-CL throughout the whole range of examined concentrations. For the latter salt, the simulated systems exhibit a drastic increase in nonideality (assuming ln fs < 0 for b = 0.5 m) upon an increase in salt concentration (i.e., the activity derivative fss is negative for b = 0.5 m and becomes more negative upon further increasing the salt concentration, hence leading to a further decrease in fs), indicating a failure of the GROMOS 54A7 and 54A8 force-field parameter sets to correctly account for the balance between ion–ion and ion–water binding propensities of sodium and chloride ions. Note in this context that the simulations performed here involve ions which are parametrized against methodology-independent hydration free energies but involve forces according to an approximate electrostatic scheme, which, owing to the predominantly negative correction terms, effectively render the ions underhydrated. It can be expected that this underhydration is more severe for cations (negative type-C1 correction) than for anions (positive type-C1 correction)[95] and possibly also more severe for monatomic ions (concentrated charge) than for oligoatomic ones (dispersed charge), which may entail an overestimation of ion pairing. Spurious clustering of ions in a 4.5 M aqueous sodium chloride solution was previously observed for the GROMOS and AMBER force fields, using ion parameters that were not calibrated with a methodologically-independent approach.[42,96]

Figure 2

Experimental salt (mean) activity coefficients fs and experimental and simulated estimates for activity derivatives fss. (a) Natural logarithm of experimental salt (mean) activity coefficients, ln(fs), displayed as a function of the natural logarithm of the salt mole fraction, ln(xs). Experimental data (circles) are from ref (79) for NA-CL and NA-ACET, ref (101) for H3C1-CL, and ref (102) for GUAN-CL. The solid lines correspond to the fit of eq A.2. (b) Experimental and simulated activity derivatives fss. The solid lines are derivatives ((∂ ln fs)/(∂ ln xs)) (eq 2) obtained from numerical centered finite differentiation of the experimental data for ln(fs) plotted in panel a. Circles denote analogous analytical derivatives at the respective salt mole fractions, obtained from fitted functions (eq A.2), using the fitting coefficients reported in Supporting Information, Table SI. Squares denote simulated results for fss (Table 2). Dashed vertical lines are a guide for the eye and indicate salt molalities bs = 0.5 m (xs = 0.0177; ln xs = −4.034), 1.0 m (xs = 0.0348; ln xs = −3.358), and 2.0 m (xs = 0.0672; ln xs = −2.700).

Cutoff artifacts turned out to significantly affect ion–water and ion–ion pair distribution functions when electrostatic interactions were handled with the Barker–Watts reaction-field scheme (data not shown). In the limit of infinite box-edge lengths, lattice-sum electrostatic interactions are Coulombic, while periodicity-induced artifacts occur with finite box-edge lengths. The authors of the present study assume that the latter artifacts are negligible for the current systems in cubic computational boxes of the size used here (about 6 nm edge length; Table 1).

In view of the nonconsideration of the properties of electrolyte solutions in the parametrization strategy of the GROMOS 54A8 force field,[22] the crude empirical relationship used to calculate heteroatomic Lennard-Jones interaction coefficients in the GROMOS force field (geometric mean combination rule),[35] and the simplicity of the employed water model (rigid three-site model without explicit polarization),[23] it is astonishing that the GROMOS 54A8 force field is capable of accounting for properties of electrolyte solutions involving the oligoatomic ions in the concentration regime up to 1.0–2.0 m. However, the present study only investigated activity coefficient derivatives. Future work will be directed toward a more thorough analysis of the thermodynamic, structural, and dynamic properties of electrolyte solutions.

Lipid Bilayers

The DPPC bilayers simulated in the present study with the GROMOS 54A8 force field (DPPC1) or with parameter sets differing solely by an increased Lennard-Jones repulsion between choline head groups and phosphate oxygen atoms (DPPC2, DPPC3) remain in the liquid crystalline-like phase, as evidenced by structural properties such as, e.g., area AL and volume VL per lipid, bilayer thickness DHH (Table 3) and electron density profiles (Figure 3a). The average area per lipid from simulation DPPC1 is at the upper threshold of experimentally observed values and equal to the reported “best” value of ref (54) (0.64 nm2), whereas increased repulsion in the bilayer headgroup moieties renders corresponding values obtained from simulations DPPC2 and DPPC3 somewhat higher (0.68 and 0.71 nm2, respectively). The average volume per lipid from simulation DPPC1 is slightly underestimated in comparison to experiment (1.222 nm3 vs 1.229–1.232 nm3), whereas corresponding values obtained from simulations DPPC2 and DPPC3 are somewhat higher (1.227 and 1.228 nm3, respectively) and thus in better agreement with experiment (Table 3). The lipid bilayer thickness deduced from the maxima of headgroup atom electron density (eq 10 and Figure 3b) across the bilayer leaflets is within the experimental range for all three DPPC simulations. Assuming an experimental value for DHH of 3.625 nm (arithmetic average between minimum and maximum values of 3.42 and 3.83 nm, respectively; Table 3), simulations DPPC1, DPPC2, and DPPC3 show deviations of 3.2, −2.3, and −5.1%, respectively. However, if the “best” value reported by ref (54) is taken as experimental reference (3.83 nm), simulation DPPC1 shows the smallest deviation (−2.3%). Considering the electron density profile due to all lipid atoms (eq 9 and Figure 3a), the distance between the maxima appears significantly smaller and underestimates the experimental value by 15.6% (DPPC1), 21.1% (DPPC2), and 24.4% (DPPC3). Tentatively approximating a measure of bilayer thickness with 2VLAL–1, the results are similar to corresponding data calculated from experimental AL and VL values. Note that the results presented above have to be taken with some caution, because (i) the simulation time scale (40 ns) might not be sufficient to allow full relaxation of lipid bilayer structural properties and exhaustive sampling of long-time scale collective events such as, e.g., undulation movements and (ii) the simulations were performed with a reaction-field electrostatic interaction function, which might cause artifacts in the sampled distribution of distances between charged groups. In the study of ref (28), some results of which are also reported in Table 3 for comparison, a production run of 120 ns length and a lattice-sum electrostatic interaction function was used to simulate a DPPC bilayer with the GROMOS 54A7 force field. Both AL and DHH agree with experimental values, whereas VL is slightly underestimated.

Table 3

Area per Lipid AL (eq 7), Volume per Lipid VL (eq 8) and Bilayer Thickness DHH Based on the Electron Density Profile Due to Lipid Head Group Atoms (eq 10 and Figure 3b) Evaluated from Simulations DPPC1, DPPC2, and DPPC3a

simulation	AL [nm2]	VL [nm3]	DHH [nm]	2VLAL–1 [nm]
DPPC1	0.64	1.222	3.74 (3.06)	3.82
DPPC2	0.68	1.227	3.54 (2.86)	3.61
DPPC3	0.71	1.228	3.44 (2.74)	3.46
C3 (54A7)b	0.63	1.226	3.57c	3.89
exp.d	0.629; 0.643	1.229; 1.232	3.42; 3.83	3.82e–3.92f
exp.g	0.64	1.232	3.83	3.85

The values reported in parentheses for DHH refer to the electron density profile due to all lipid atoms (eq 9 and Figure 3a). For comparison, corresponding results obtained by ref (28) with the GROMOS 54A7 force field are also reported. Experimental data are provided in the form of extreme values found after an investigation of various data sources (see also ref (28)) along with “best” values reported in the review of ref (54).

Results of simulation C3 reported by ref (28) using the GROMOS 54A7 force-field parameter set and a production run of 120 ns.

Ref (28) seems to define the position of the head groups from the total hydrated lipid bilayer electron density profile.

Extrema of experimental values reported by refs (52−55), refs (53−55), and refs (53−55,103) for AL, VL, and DHH, respectively.

Value obtained based on AL = 0.643 nm2 and VL = 1.229 nm3.

Value obtained based on AL = 0.629 nm2 and VL = 1.232 nm3.

“Best” values reported in the review of ref (54).

Figure 3

Electron density profile for the DPPC bilayer in simulations DPPC1, DPPC2, and DPPC3. The former corresponds to the GROMOS 54A8 force-field parameter set, whereas the latter two simulations involve increased Lennard-Jones repulsion between the choline head groups and the phosphate oxygen atoms. The dashed lines indicate extreme values for experimentally observed bilayer thicknesses (minimum of[103] 3.42 nm and maximum of[54] 3.83 nm, respectively; see also ref (28) for a compilation of experimental data sources). (a) Electron density profile ηe(z) due to all lipid atoms (eq 9). (b) Electron density profile ηe(H)(z) due to lipid headgroup atoms (eq 10).

The Lennard-Jones repulsion coefficient C12, for interactions between the united-atom methyl groups of the choline head groups and phosphate oxygen atoms evaluates to 6.93 × 10–6, 1.48 × 10–5, 2.47 × 10–5, and 1.58 × 10–5 kJ·mol–1·nm12 in simulations DPPC1, DPPC2, and DPPC3 and the GROMOS 54A7 force field, respectively, so that actually the force-field parameter set employed in simulation DPPC2 is most similar to the GROMOS 54A7 force field. However, the results obtained in the present study (Table 3) give better agreement with experimental values in the case of DPPC1 concerning AL, DHH and 2VLAL–1. This difference is most likely caused by the usage of a different electrostatic interaction function. It has been pointed out before that AL is very sensitive to simulation methodology, in particular, to the treatment of long-range electrostatic interactions.[97] Notably, larger AL values were found before with a reaction-field electrostatic interaction function than with a lattice-sum electrostatic interaction function.[97] Since it is difficult to assess the relative merits of different electrostatic interaction functions in the absence of an “ideal” reference characterized by macroscropic nonperiodic extent and Coulombic electrostatic interactions, and since the parametrization of the GROMOS force field relies on a truncated electrostatic interaction function with a Barker–Watts reaction-field correction, simulation DPPC1 is still considered to provide the most reasonable area per lipid in the present context.

Time evolutions of AL and VL during the production runs are shown in Supporting Information Figures S1 and S2. While simulations DPPC2 and DPPC3 present a sharp initial increase in AL, simulation DPPC1 shows a rather stable area per lipid around 0.64 nm2. The volume per lipid VL exhibits less pronounced fluctuations than AL. It has been suggested before that VL might therefore be a more suitable property for the purpose of force-field calibration or/and validation.[28] However, its calculation involves an estimate for the effective volume of a water molecule (eq 8), evaluated from the simulation of a pure water system. The authors of the present study think that it might not be appropriate to use such an estimate for water in the presence of a lipid bilayer, because solvent electrostriction effects in the vicinity of the charged head groups will not be properly accounted for.

Proteins

Secondary Structure Stability

As evidenced by quantitative structural characteristics of secondary structure elements, such as, e.g., the number of backbone hydrogen bonds Nbb, the rmsd of the backbone from a representative reference structure rmsdbb, or the fraction of residues present in given secondary structure elements, the proteins simulated in the present study remain structurally intact (Table 4). This is not surprising because (i) the suitability of the GROMOS 54A7 force field for protein simulation has been previously validated[21] and (ii) the sole difference between the GROMOS 54A7 and 54A8 force fields is in the description of charged amino acid side chains, and it is known that intramolecular salt bridges influence protein stability only marginally.[43−46] As can be seen from Table 4, the average backbone rmsd values rmsdbb of all proteins are within 0.36 nm of the initial structure. PROTG and COLDS exhibit exceptional stability, with rmsdbb values of 0.13 and 0.17 nm, respectively. The highest value is found for SACNaCl54A8 (0.36 nm) and can be drawn back to a movement of the last ca. seven C-terminal residues (see below). For FOX, the reported rmsdbb value refers to residues 9–82 only. This protein has extremely flexible terminal tails, and the rmsdbb calculated for the whole protein (residues 1–88) amounts to 0.54 nm. The time evolutions of rmsdbb values are provided in Supporting Information Figure S3 and also illustrate that deviations from the initial structure in FOX can mainly be ascribed to structural fluctuations occurring at the termini.

On average, the initial number of backbone hydrogen bonds is maintained up to 86.7% in simulations HEWL, FOX, CM, GCN, PROTG, and COLDS, respectively, and 78.4% in simulations SAC054A7, SAC054A8, SACCl54A7, SACCl54A8, SACNaCl54A7, and SACNaCl54A8, respectively. These values are relatively high (>75%), i.e., indicative of an overall secondary structure close to the initial one, with the exception of SACNaCl54A8 where the C-terminal tail formed by residues 60–66 (RAEREKK) moves at the beginning of the simulation. This is reflected by an increased root-mean-square fluctuation (rmsf) of the corresponding Cα atoms (the average Cα rmsf values of residues 1–59 and 60–66 are 0.19 and 0.40 nm, respectively), an associated loss of α-helical structure (Figure 5), and the fact that the average backbone rmsd obtained upon the omission of the seven terminal residues from the calculation is lower than the value of 0.36 nm reported in Table 4, and more similar to the values found for the other SAC simulations, namely 0.31 nm. Visual inspection of the coordinate trajectory suggests that the tail movement is accidental (or water-driven; see, however, below), because it could not be correlated with any other obvious event such as, e.g., counterion binding. The possibly most crucial change between the GROMOS 54A7 and 54A8 force fields affects the carboxylate side chains, which are more hydrophilic by about 40 kJ·mol–1 in the latter parameter set. Therefore, one might feel tempted to ascribe an increased flexibility of the tail to the hydration properties of the glutamate side chains E62 and E64 and the carboxylate terminus of K66. This is, however, unlikely, because the α-helical tail structure is e.g. maintained throughout simulation SACCl–54A8 and is affected by similar stability problems toward the end of simulations SAC054A7 and SACCl–54A7 (Figure 5).

Figure 5

Residues occurring in secondary structure elements α-helix (black), 310-helix (green), π-helix (blue), or β-sheet (red) during the simulations of SAC054A7, SAC054A8, SACCl54A7, SACCl54A8, SACNaCl54A7, and SACNaCl54A8. Secondary structure elements of the initial (minimized) structure are indicated in the short stretch at the left-hand side of each graph.

Besides monitoring backbone hydrogen bonds, other indications of stable secondary structure are the average fractions of protein residues in given secondary structure elements. Overall, these values are also close to those characterizing the initial structure. In combination with visual inspection of the coordinate trajectories or/and the time evolution of secondary structure elements (Figures 4 and 5), this suggests the absence of a severe spurious structural bias in the underlying force field. The most striking secondary structural changes are an increase in α-helix content of HEWL and CM; a decrease in α-helix content of GCN and SAC; a decrease in 310-helix content of HEWL, CM, GCN, and SAC; an increase in π-helix content of GCN; an increase in β-sheet content of COLDS; and a decrease in β-sheet content of SAC. Clearly, in HEWL and CM, there is transformation between 310- and α-helices, and in GCN, there is a transformation between a π- and an α-helix. The 310-helices formed by residues 120–123 of HEWL and 2–3 of CM convert to α-helices after about 1.7–2.2 and 0.6 ns, respectively, and largely retain this conformation during the remainder of the simulation time. It was suggested before[21] that the GROMOS 54A7 force field marginally overstabilizes α-helices in comparison to 310-helices, and the same finding appears to hold for the GROMOS 54A8 force field. On the contrary, the transformation to a π-helix in the GCN peptide after about 20 ns seems to be more transient since the α-helical conformation is recovered after about 45 ns. For FOX, one can recognize a loss of β-sheet (residues 71–73) after about 7.4 ns. This is not severe because the compatibility with experimental NOE values is still ensured (see below).

Figure 4

Residues occurring in secondary structure elements α-helix (black), 310-helix (green), π-helix (blue), or β-sheet (red) during the simulations of HEWL, FOX, CM, GCN, PROTG, and COLDS. Secondary structure elements of the initial (minimized) structure are indicated in the short stretch at the left-hand side of each graph.

The SAC protein seems to be somewhat less stable than the other proteins, which can, to a certain extent, be drawn back to the flexibility of the last ca. seven C-terminal residues, which is also qualitatively reflected in the reported “B-factors” of the NMR structure. Note also that the SAC protein is native to a hyperthermophilic organism (Table 1). At elevated temperatures, the solvation properties of water are considerably different from those at 300 K, i.e., due to a decrease in relative dielectric permittivity, water becomes less solvating. This enhances protein stability at high temperatures (enhanced salt-bridge formation). A limited loss of secondary structure elements at 300 K as encountered in the present study might therefore not be observed in simulations under the physiological conditions of the hyperthermophile species Sulfolobus acidocaldarius. However, no such simulations were undertaken because the temperature dependence of the relative dielectric permittivity of the SPC water model does not quantitatively reproduce the experimental behavior (data not shown). The structural flexibility observed in MD simulations might merely be in keeping with the common notion that structures deposited in the PDB are not necessarily representative of the entire configurational ensemble sampled by the protein in the course of time.

In the SAC protein (Figure 5), a loss of α-helical structure at the C-terminus is found for SAC054A7, SACCl54A7, SAC054A8, and SACNaCl54A8. It is likely due to the flexibility of the termini, which is not accounted for by the single structure deposited in the PDB. Moreover, all SAC simulations present a transient formation/deformation of a 310-helix (residues 47–49) present in the NMR structure. In simulation SACCl54A7, another 310-helix is similarly destabilized (residues 17–19). No clear-cut differences can be observed depending on the counterion concentration. Overall, the SAC protein seems to have more pronounced flexibility in its secondary structure than the other proteins examined in the present study.

Experimental NMR data in the form of NOE values and 3J-coupling constants is available for HEWL, FOX (NOE values only), and GCN. For all three proteins, the upper thresholds derived from the NOE values are satisfied to a good extent during the simulations (Figure 6a, c, d), corresponding violations being similar to those observed previously in simulations using the GROMOS 54A7 force field[21] (note the erroneous reporting of violations for GCN in Figure 11 of ref (21)). Furthermore, it can be seen that the loss of β-sheet in the FOX protein after about 7.4 ns does not affect the quality of NOE upper bound fulfillment. Nine (first 7.4 ns) or 18 (whole simulation) of 79 NOE values involving at least one of residues 71, 72, and 73 show violations, with average and maximum values of 0.011 and 0.25 nm, respectively, during the first 7.4 ns and 0.018 and 0.23 nm, respectively, during the whole simulation. Violations of experimental 3J-coupling constants between backbone hydrogen atoms are of similar magnitude to those observed previously in simulations using the GROMOS 54A7 force field[21] (note the erroneous reporting of violations for HEWL in Figure 5 of ref (21)). In general, 3J-coupling constants are reproduced less well than NOE values (Figure 6b, e), which is mostly due to (i) the empirical nature of the relationship between simulated dihedral angles and corresponding deduced 3J-coupling constants (eq 12), giving rise to uncertainties of up to 1 Hz in the calculated values; (ii) insufficient configurational sampling in the present simulations in comparison to the time scales experimental 3J-coupling constants are averaged over; (iii) inaccurate experimental data in the case of GCN.[89] In particular, as attested by the investigation of backbone hydrogen bonds and secondary structure stability (Table 4 and Figure 4), proteins HEWL and GCN appear stable despite the indication of severe violations of experimental data by the backbone hydrogen 3J-coupling constants.

Figure 6

Comparison of the simulations of HEWL (a, b), FOX (c), and GCN (d, e) with experimental NMR data in terms of NOE values (a, c, d) and 3J-coupling constants (b, e; HEWL and GCN only). Bar diagrams depict the fractions of violations in bins of 0.05 nm (NOE values) or 0.5 Hz (3J-coupling constants), where corresponding bars are drawn at the bin center locations. The fractions rely on a total of 1630, 1518, and 179 NOE values for HEWL, FOX, and GCN, respectively, as well as 95 and 15 3J-coupling constants for HEWL and GCN, respectively. The lines represent cumulative sums of the fractions of violations. For FOX, additional cumulative sums obtained from the first 7.4 ns of the simulation are shown in green color. Note that for NOE values, violations refer to positive deviations, whereas for 3J-coupling constants, violations refer to absolute deviations from the experimental data. For the latter, deviations with a magnitude of less than 0.1 Hz were defined as zero deviations (three occurrences).

Salt-Bridge Formation and Counterion Binding Properties

As it is single-ion hydration free energies of charged amino acid side chain analogs which distinguish the GROMOS 54A8 force field from its predecessor 54A7, the thermodynamic equilibrium between salt-bridge formation and hydration of counterions as well as of charged functional groups at the protein surface might be shifted. Possible differences would be reflected in the configurational sampling of distances between charged groups and between charged groups and the solvent, as well as of the orientations of these moieties with respect to each other. However, counterion distributions around (bio)molecules are often difficult to converge due to relatively slow diffusion of the ions. This is why explicit addition of counterions is sometimes omitted, as, e.g., done in a previous simulation study of the SAC protein.[98] The authors of the present study deem the performed 20 ns simulations sufficiently long to assess possible differences in configurational sampling of charged groups between the GROMOS 54A7 and 54A8 force-field descriptions.

Configurational sampling of the ionic and solvent groups of interest can be characterized by probability distribution functions or/and the averages and fluctuations of corresponding time series. Here, distance probability distributions between charged groups and between charged groups and the solvent were chosen to detect possible consequences of the alternate representation of charged amino acid side chains. Figure 7 depicts radial distribution functions g(r) (eq 13) of the number density of water oxygen atoms, sodium ions, and chloride ions around arginine, lysine, aspartate, and glutamate residues of the SAC protein in simulations SACCl54A7, SACCl54A8, SACNaCl54A7, and SACNaCl54A8. It can be seen that using the GROMOS 54A8 force field in comparison to the GROMOS 54A7 force field (i) the water number density is enhanced around lysine, aspartate, and glutamate residues, water molecules approach the charged functional groups of arginine and lysine residues more closely, and the first hydration shell of the terminal guanidinium group in arginine residues is tighter, i.e., the corresponding peak is less broad; (ii) the number density of sodium ions is reduced around arginine residues and enhanced around aspartate and glutamate residues, whereas it is very similar around lysine residues; (iii) the number density of chloride ions is enhanced around arginine and lysine residues and decreased around aspartate and glutamate residues, and chloride ions approach the charged functional groups of arginine and lysine residues more closely. These effects are to be expected based on the stronger hydration of guanidinium and acetate ions in the GROMOS 54A8 force field (due to charge enhancement). The situation is somewhat more complex in the case of the lysine side chain, the small-molecule analog of which is the methylammonium ion in the parametrization procedure. The hydration free energy of this ion is (approximately) the same in the GROMOS 54A7 and 54A8 force-field parameter sets. Nevertheless, the force-field representation of this ion is not the same in both parameter sets. While the distribution of partial atom charges is rather similar, the methyl group Lennard-Jones radius was increased whereas the nitrogen Lennard-Jones radius was decreased upon introduction of the GROMOS 54A8 force-field version.[22] The lysine side chain in the GROMOS 54A8 force field did not adopt a larger Lennard-Jones radius for the methylene group attached to the side-chain nitrogen atom. Therefore, the lysine side chain might in reality be somewhat overhydrated. This explains the above findings of increased water and chloride ion number density around this residue.

Figure 7

Radial distributions g(r) (eq 13) of the number density of water oxygen atoms (OW), sodium ions (NA), and chloride ions (CL) around arginine (black), lysine (red) as well as aspartate and glutamate (green) residues of the SAC protein in simulations involving solely a neutralizing amount of chloride counterions (SACCl–) or a neutralizing amount of chloride counterions in combination with an equimolar sodium chloride solution (SACNaCl). The distance measurements are based on the nitrogen atoms of the alkylguanidinium group of arginine residues, the nitrogen atoms of the alkylammonium group of lysine residues, and the carboxylate oxygen atoms of aspartate and glutamate residues. A bin width of 0.01 nm was used in the calculation of g(r). Solid and dashed lines refer to simulations with the GROMOS 54A7 (SACCl54A7, SACNaCl54A7) and 54A8 (SACCl54A8, SACNaCl54A8) force field, respectively.

Figure 8 depicts normalized (eq 14) probability distribution functions Pαβ(r) of the (minimum-image) distance of all pairs of like-charged positive, like-charged negative, and oppositely-charged protein functional groups of the SAC protein in simulations SAC054A7, SAC054A8, SACCl54A7, SACCl54A8, SACNaCl54A7, and SACNaCl54A8. The curves obtained from the two force-field parameter sets are very similar. Considering first peak positions and heights, only marginal differences can be observed, i.e., using the GROMOS 54A8 force field in comparison to the GROMOS 54A7 force field, it appears that the first peak of the distance probability distribution (i) between like-charged positive protein functional groups is higher in simulations without counterions or solely a neutralizing amount of chloride counterions, whereas peak positions are unchanged; (ii) between like-charged negative protein functional groups is somewhat higher and followed by an immediate second peak in the simulation without counterions, is smaller in the simulation with solely a neutralizing amount of chloride counterions and is smaller but occurs at a considerably shorter distance in the simulation containing a neutralizing sodium chloride solution; (iii) between oppositely-charged protein functional groups is not affected. Note that, overall, the shape of the probability distributions involving like-charged positive and oppositely-charged protein functional groups are extremely similar, whereas the probability distributions involving like-charged negative protein functional groups present more variation between the two different force fields, especially in the simulations containing a neutralizing sodium chloride solution (SACNaCl54A7, SACNaCl54A8). The SAC protein carries a net charge of +6e (Table 4), and therefore, presumably, probability distribution features concerning the separation of positive charges and the formation of salt-bridged ion pairs converge quicker than features concerning the separation of negative charges. Due to possible convergence issues, the latter should therefore be taken with some caution. The equivalence of the distance probability distributions involving oppositely-charged protein functional groups, in particular, the occurrence of an equally high salt-bridged ion-pair peak at a distance of 0.56 nm, is not surprising because potential of mean force calculations as a function of interionic separation of the corresponding charged side chain analogs in water were already found to be rather similar.[22]

Figure 8

Probability distributions Pαβ(r) (section II.3.2) of the (minimum-image) distance of all pairs of like-charged positive (α = +, β = +), like-charged negative (α = −, β = −), and oppositely-charged (α = +, β = −) protein functional groups of the SAC protein in simulations involving no counterions (SAC0), solely a neutralizing amount of chloride counterions (SACCl–), or a neutralizing amount of chloride counterions in combination with an equimolar sodium chloride solution (SACNaCl). The functions are normalized to integrate to the number of distinguishable pairs (eq 14). The distance measurements are based on the alkylguanidinium carbon atoms of arginine residues, the alkylammonium carbon atoms of lysine residues, and the carboxylate carbon atoms of aspartate and glutamate residues. A bin width of 0.01 nm was used in the calculation of Pαβ(r). Black and red lines refer to simulations with the GROMOS 54A7 (SAC054A7, SACCl54A7, SACNaCl54A7) and 54A8 force field (SAC054A8, SACCl54A8, SACNaCl54A8), respectively. The dashed vertical lines indicate an interionic distance of 1.4 nm, i.e., the cutoff distance employed during configurational sampling with the Barker–Watts reaction-field scheme (section II.2).

Finally, it should be pointed out that distance probability distributions involving charged species will be affected by the use of an approximate electrostatics scheme. A cutoff-truncated effective electrostatic interaction function with a Barker–Watts reaction-field correction was previously shown to provoke the appearance of a spurious maximum just below and of a spurious minimum just above the cutoff distance in the potential of mean force between oppositely-charged ion pairs in a homogeneous dielectric of a relative permittivity corresponding to that of water.[99,100] The same phenomenon is observed in Figure 8, where a minimum and a maximum in the distance probability distribution between oppositely-charged protein functional groups are observed just below and above the cutoff distance of 1.4 nm, respectively. Similar artifacts can be found in the potential of mean force between like-charged ion pairs.[99,100] Most importantly, the cutoff-induced artifacts turned out to be of significant magnitude in comparison to the thermal energy,[99,100] and it should thus be kept in mind that the distance probability distributions presented in Figure 8 are likely to involve artifacts in side-chain configurational sampling due to the usage of approximate electrostatic interactions.

Conclusion

The present study investigated the ability of the GROMOS 54A8 force field to accurately model the structural properties of lipid bilayers, proteins, and electrolyte solutions. Its aim was 4-fold (points i–iv in section I), and the conclusions reached in this study can be summarized as follows:

(i) Concerning aqueous electrolyte solutions involving ionic species parametrized against methodology-independent hydration free energies, a realistic description of the thermodynamic equilibrium between ion–ion and ion–water pairing propensities in the GROMOS 54A8 force field is not a priori guaranteed due to electrostatic artifacts in configurational sampling, the crude empirical relationship used to calculate heteroatomic Lennard-Jones interaction coefficients in the GROMOS force field, and the simplicity of the employed water model. Here, it was assumed that the present simulations of electrolyte solutions with lattice-sum electrostatic interactions are free of periodicity-induced artifacts in cubic computational boxes of about 6 nm edge length. Salts involving oligoatomic species (NA-ACET, GUAN-CL, H3C1-CL) are found to reproduce experimental salt activity derivatives for concentrations up to 1.0 m very well, and good agreement between simulated and experimental data is also found for NA-ACET and H3C1-CL at 2.0 m concentration. However, not even qualitative agreement is found for NA-CL throughout the whole range of examined concentrations, indicating a failure of the GROMOS 54A7 and 54A8 force-field parameter sets to correctly account for the balance between ion–ion and ion–water binding propensities of sodium and chloride ions.

(ii) Concerning lipid bilayers, the GROMOS 54A8 force field differs from its predecessor 54A7 by a modified Lennard-Jones repulsion between the choline head groups and the phosphate oxygen atoms. It is found that the GROMOS 54A8 force field reproduces the liquid crystalline-like phase of a hydrated DPPC bilayer at a pressure of 1 bar and a temperature of 323 K. The area per lipid is in agreement with experimental data, whereas other structural properties (volume per lipid, bilayer thickness) appear slightly underestimated. Considering area per lipid estimates from MD simulation, the strong dependence on simulation methodology, in particular, the treatment of long-range electrostatic interactions, should be kept in mind. The present study used truncated electrostatic interactions with a Barker–Watts reaction-field correction for the simulation of lipid bilayers, and it is likely that a slightly lower area per lipid value will be obtained with a lattice-sum electrostatic interaction function.

(iii) Concerning proteins, the GROMOS 54A8 force field differs from its predecessor 54A7 in the description of charged amino acid side chains. A range of different proteins was simulated at pH 7, i.e., with all arginine, lysine, aspartate, and glutamate residues carrying net charges. The secondary structure of these proteins is largely maintained and found compatible with experimental NMR data (NOE values and 3J-coupling constants; where available), which is not surprising because a putatively different thermodynamic equilibrium between salt-bridge formation and hydration of charged functional groups at the protein surface introduced by the new description of charged amino acid side chains is expected to influence protein stability only marginally.[43−46] It should be pointed out that, as was observed for the GROMOS 54A7 force field,[21] the GROMOS 54A8 force field appears to slightly overstabilize α-helices in comparison to 310-helices.

(iv) Concerning proteins in aqueous electrolyte solutions, it appears also important to investigate whether usage of the GROMOS 54A8 force field leads to an altered counterion binding behavior and, as a consequence thereof, a possible impact on salt-bridge formation between charged side chains at the protein surface. On the basis of the law of “matching water affinities”,[40,47,48] it cannot a priori be ruled out that the counterion distribution around charged protein residues will be affected by the different hydration properties of these residues. Indeed, slight differences in the radial distribution functions of water molecules, sodium and chloride counterions around protein charged functional groups are detected, which can be explained based on the different hydration properties of the corresponding ionic side chain analogs. On average, the side chains of arginine, lysine, aspartate, and glutamate residues appear slightly more hydrated and present a slight excess of oppositely-charged solution components in their vicinity. Salt-bridge formation properties between charged residues at the protein surface, as assessed by probability distributions of interionic distances, are largely equivalent in both force fields. Considering configurational sampling of charged moieties in MD simulation, the strong dependence on the treatment of long-range electrostatic interactions should be kept in mind.

Overall, the authors of the present study think that the 54A8 parameter set of the GROMOS force field can be applied in (bio)molecular simulations of electrolyte solutions, lipids, and proteins, as well as mixtures thereof. Comparison in the context of protein and lipid simulations with its predecessor 54A7 suggests comparable performance. Since the charged functional groups in the GROMOS 54A8 parameter set are derived from a rigorous comparison to single-ion hydration free energies, this parameter set is considered an important step in the continuous development of the GROMOS force field.