Evaluating the strength of salt bridges: a comparison of current biomolecular force fields.

Recent advances in computer hardware and software have made rigorous evaluation of current biomolecular force fields using microsecond-scale simulations possible. Force fields differ in their treatment of electrostatic interactions, including the formation of salt bridges in proteins. Here we conducted an extensive evaluation of salt bridge interactions in the latest AMBER, CHARMM, and OPLS force fields, using microsecond-scale molecular dynamics simulations of amino acid analogues in explicit solvent. We focused on salt bridges between three different pairs of oppositely charged amino acids: Arg/Asp, Lys/Asp, and His(+)/Asp. Our results reveal considerable variability in the predicted KA values of the salt bridges for these force fields, as well as differences from experimental data: almost all of the force fields overestimate the strengths of the salt bridges. When amino acids are represented by side-chain analogues, the AMBER ff03 force field overestimates the KA values the least, while for complete amino acids, the AMBER ff13α force field yields the lowest KA value, most likely caused by an altered balance of side-chain/side-chain and side-chain/backbone contacts. These findings confirm the notion that the implicit incorporation of solvent polarization improves the accuracy of modeling salt bridge interactions.

Introduction

Recent advances in computer hardware and software have greatly extended the time scales that can be covered by biomolecular simulations. These longer time scales (beyond nanoseconds) are essential for the rigorous evaluation of current biomolecular force fields. One important characteristic of these force fields is the ability to accurately model the formation of salt bridges, or pairs of amino acids whose oppositely charged side-chains are within hydrogen-bonding distance in proteins.[1] However, it has long been suspected that the forces between oppositely charged amino acids are overly attractive in molecular dynamics (MD) simulations with current biomolecular force fields, and there have been a number of efforts to reduce this artifact in the improvement of various force fields.[2−4] Previous theoretical studies have analyzed the contribution of salt bridges to protein or protein–protein complex stability, using both implicit[5−9] and explicit modeling of solvation.[10,11] Others have studied salt bridges using amino acid analogues,[12−20] often employing biasing techniques in the simulations.[13−15] More recently, a comprehensive comparison of force field/water model combinations was conducted for salt bridge interactions between the amino and carboxyl groups of zwitterionic amino acids, using extensive simulations in explicit solvent on the microsecond time scale.[21]

Here, we evaluated six biomolecular force fields for their ability to accurately model the strengths of salt bridges between the side-chains of oppositely charged amino acids by unbiased, microsecond-scale MD simulations in explicit solvent. In particular, we directly compared current AMBER, CHARMM, and OPLS force fields in simulations of association between the side-chain analogues of three different pairs of amino acids, Arg/Asp, Lys/Asp, and His(+)/Asp. We further tested one of the pairs, Arg/Asp, by simulating association of blocked amino acid dipeptides. In addition, we evaluated the influence of the solvent model on the strengths of the salt bridges by simulating the side-chain analogue pairs using a selection of different force field/water model combinations. To our knowledge, our microsecond-scale simulations provide the most extensive sampling of salt bridge formation to date, yielding thousands of association/dissociation events, permitting quantitative comparisons, both between the force fields and with experiment. Our results reveal considerable variability among the current force fields in terms of the resulting strengths of salt bridge interactions, as well as differences from experimental data.

Methods

Preparation of Starting Models

We modeled the formation of salt bridges between the following pairs of oppositely charged amino acids using side-chain analogues: Arg/Asp (guanidinium cation/acetate anion), Lys/Asp (butylammonium cation/acetate anion), and His(+)/Asp (imidazolium cation/acetate anion). Our systems were constructed to be consistent with the experimental conditions under which the equilibrium association constants (KA) of guanidinium acetate and butylammonium acetate have been measured,[22] i.e., using the same concentrations (0.9 M guanidinium and 0.02 M acetate, which corresponds to 100 molecules of guanidinium and two molecules of acetate in the presence of ∼18 000 explicit water molecules). To ensure a net charge of zero, we included 98 chloride ions (the same counterion that is present in the experiments). The same concentrations of the cation, anion, and chloride ions were also used for the model systems consisting of butylammonium/acetate and imidazolium/acetate. Starting models for these simulations were constructed using the Packmol software package,[23] immersing the appropriate number of side-chain analogues in periodic, cubic boxes of explicit solvent. For the Arg/Asp salt bridge, we also used blocked amino acid dipeptides (acetyl—arginine—N-methyl and acetyl—aspartate—N-methyl) to model salt bridge formation. Only a single copy of each blocked dipeptide was included, corresponding to a concentration of 0.012 M for each salt-bridging partner, with a distance of 10 Å between the amino acids. All force field parameters of the side-chain analogues were based on those of the complete amino acids. For the chloride ions, parameters derived specifically for the water model were used when available; otherwise, parameters derived for a similar water model were used.[24−26] Nonbonded parameters of the side-chain analogues, along with those used to model chloride ions and blocked amino acid dipeptides, are provided in Table S1 of the Supporting Information.

To alleviate any unfavorable interactions, each model was subjected to energy minimization followed by a two-stage equilibration with harmonic position restraints on all heavy atoms of the side-chain analogues (force constant of 10 kcal mol–1 Å–2) using the Desmond 3.0.1.0 software package.[27] In the first stage, the energy-minimized system was equilibrated for 20 ps at constant temperature (25 °C) using a weak Langevin thermostat (frictional constant of 1 ps–1). During the second stage, the system was equilibrated for 1 ns at constant temperature (25 °C) and pressure (1 atm) using the Martyna–Tobias–Klein thermostat and barostat[28] (coupling time constants of 1.0 and 2.0 ps, respectively). To enable a 2 fs time step, bonds to hydrogen were constrained to their equilibrium values using the M-SHAKE algorithm.[29] A short-range nonbonded cutoff of 10.0 Å was used, and long-range electrostatics were calculated using the particle mesh Ewald (PME) method.[30] The frame from the second half of the NPT equilibration with volume closest to the average was used to start the production simulation.

Simulation Details

To obtain extensive sampling of salt bridge association (and dissociation) events, 1-μs MD simulations were performed for each side-chain analogue system; 10-μs simulations were performed for the blocked arginine and aspartate dipeptide systems. All simulations were carried out in the NVT ensemble using a 64-node Anton special-purpose supercomputer, which is able to run MD simulations roughly 2 orders of magnitude faster than conventional hardware[31] (altogether, the simulations required a total of 40 machine-days). The temperature was maintained at 25 °C using the Nosé–Hoover thermostat with a weak coupling constant of 0.5 ps.[32] van der Waals and short-range electrostatic interactions were truncated at 10.0 Å; long-range electrostatic interactions were calculated using the Gaussian split Ewald method.[33] To enable a 2.5 fs time step, bonds to hydrogen were constrained to their equilibrium lengths using the M-SHAKE algorithm.[29] Conformations were saved every picosecond for analysis.

Calculation of Equilibrium Association Constants

Equilibrium association constants (KA) were calculated from the populations of the bound and unbound states of the oppositely charged side-chain analogues. For example, the KA for association between guanidinium and acetate was calculated using the following:where Pbound is the population of the bound state, Punbound guanidinium and Punbound acetate are the populations of unbound guanidinium and acetate, respectively, and C0 is the reference concentration of guanidinium (i.e., 0.9 M). In addition to species in which a single acetate molecule is bound to a single cation molecule, forming a 1:1 complex (e.g., the guanidinium/acetate complex), species in which acetate is bound to two cation molecules, forming a 1:2 complex (e.g., the diguanidinium/acetate complex) were observed. KA values for the latter are included in Table S2 of the Supporting Information; the results discussed below focus on formation of the major complex, which is the 1:1 complex. Standard errors in the KA values were calculated using a block averaging method.[34]

For each side-chain analogue system, the unbound and bound states were defined using the potential of mean force (PMF) as a function of the minimum distance between the nitrogen and oxygen atoms of the positively and negatively charged analogues, respectively (minimum N–O distance; see Figure 1). In particular, the point of inflection between the bound state free energy minimum (∼2.5–3 Å) and the desolvation barrier (∼3–3.5 Å) was used as the bound state cutoff, while 4.5 Å was used as the unbound state cutoff. If the minimum N–O distance between an analogue pair dropped below the bound state cutoff they were classified as bound until they crossed the unbound state cutoff, and vice versa. For simulations of the blocked arginine and aspartate dipeptides, the same definitions of the unbound and bound states were used as for the guanidinium/acetate system.

Figure 1

Potentials of mean force (PMF) between three different pairs of oppositely charged side-chain analogues using six biomolecular force fields with the TIP4P-Ew explicit water model.

Calculation of the Solvent Dielectric Constant

The dielectric constant of water in each simulation, εwater, was calculated using the following equation:where Mwater is the net dipole moment of water, Vwater is the volume occupied by water, T is the temperature of the system, kB is the Boltzmann constant, and ε0 is the permittivity of free space. The net dipole moment of water was calculated using the following:where q is the atomic charge of each water site, and r is its position vector.[35,36] The dielectric constant of water, rather than that of the complete system, was used since it is impossible to calculate the contributions of molecules with a net charge to the system dipole moment from simulations with periodic boundary conditions.[36] The appropriate volume was thus the volume of the water molecules present in the system. For each water model used, a pure water system of the same total volume as that of the side-chain analogue systems was equilibrated using the same protocol, and the molecular volume of water calculated. For each analogue system, the number of water molecules present was multiplied by the molecular volume to calculate the approximate volume of water present in the system. Standard errors in the εwater values were calculated using a block averaging method.[34]

Results and Discussion

We compared the following six current biomolecular force fields in terms of their ability to model salt bridge interactions: AMBER ff99SB-ILDN,[37] AMBER ff03,[2] AMBER ff13α,[4] CHARMM27,[38] CHARMM22*,[3] and OPLS_2005.[39] In particular, we simulated association (and dissociation) of salt bridges between the following three pairs of oppositely charged amino acids: Arg/Asp, Lys/Asp, and His(+)/Asp. We focused primarily on simulating side-chain analogues (i.e., guanidinium, butylammonium, and imidazolium cations for arginine, lysine, and histidine, respectively, and acetate anion for aspartate) since these analogues are the minimal systems for studying the formation of salt bridges; in addition, equilibrium association constant (KA) values for such systems have been experimentally measured, providing an excellent opportunity to validate the simulations. While blocked amino acid dipeptides (i.e., acetyl—amino acid—N-methyl) might be regarded as being more representative of the protein environment, no experimental KA values for the association of oppositely charged amino acid dipeptides are available. Nonetheless, we evaluated the force fields in simulating such systems, focusing on just one of the three salt bridges, Arg/Asp. Finally, in addition to the above simulations, in which each biomolecular force field was paired with the TIP4P-Ew explicit water model, which reproduces the liquid properties of water at the temperatures and pressures relevant to biology,[40] we also evaluated the influence of the water model on the strength of the salt bridges by testing a selection of force field/water model combinations for all three pairs of side-chain analogues. For each force field a selection of water models drawn from TIP3P,[41] mTIP3P,[42] TIP4P,[41] TIP4P/2005,[43] and SPC/E[44] were tested, including the water model with which each force field was originally derived.

Association Constants of Side-Chain Analogues

To validate our simulations of association between oppositely charged side-chain analogues, we computed KA values and compared these to those measured by experiments. Experimental KA values have been measured for guanidinium/acetate and butylammonium/acetate association by monitoring changes in the pKa of acetate in the presence or absence of either the guanidinium or butylammonium cation.[22] Our microsecond-long simulations yielded thousands of independent binding events, permitting the extraction of extremely precise KA values, with the mean lifetimes of the bound state ranging from ∼10–300 ps and the mean lifetimes of the unbound state ranging from ∼20–120 ps (Table S2, Supporting Information).

In general, the KA values computed from our side-chain analogue simulations are overestimated in comparison to experimentally measured values, with the AMBER ff03 force field overestimating the strengths of the salt bridges to the least extent and the OPLS_2005 force field to the greatest extent, when using the same water model (Table 1). The computed KA values vary considerably among the force fields. For example, when using the TIP4P-Ew water model, the computed KA values for the three types of salt bridges vary by as much as ∼4-fold, ∼3-fold, and ∼4-fold for the associations of guanidinium, butylammonium, and imidazolium with acetate, respectively, which amounts to ∼1.4-fold, ∼1.8-fold, and ∼1.9-fold differences in the probabilities of binding (Pbound) (see Figure 2). We note that our definition of the bound state is very conservative and that the use of less conservative definitions (e.g., use of the desolvation barrier as a cutoff) yields even stronger association constants, without affecting our overall conclusions.

Table 1

Association Constants (KA) and Probabilities of Binding (Pbound) for Three Different Pairs of Oppositely Charged Side-Chain Analogues Using Six Biomolecular Force Fields and Six Explicit Water Modelsa

		guanidinium/acetate		butylammonium/acetate		imidazolium/acetate
force field	water model	Pboundb	KA (M–1)	Pboundb	KA (M–1)	Pboundb	KA (M–1)
AMBER ff99SB-ILDN	TIP4P-Ew	0.57	2.23 ± 0.03	0.32	0.53 ± 0.01	0.36	0.65 ± 0.01
AMBER ff03	TIP4P-Ew	0.45	1.12 ± 0.01	0.28	0.45 ± 0.00	0.32	0.54 ± 0.00
AMBER ff13α	TIP4P-Ew	0.54	2.28 ± 0.03	0.48	1.20 ± 0.01	0.40	0.79 ± 0.01
CHARMM27	TIP4P-Ew	0.61	4.06 ± 0.10	0.45	0.98 ± 0.01	0.35	0.63 ± 0.01
CHARMM22*	TIP4P-Ew	0.48	1.31 ± 0.02	0.39	0.75 ± 0.01	0.29	0.47 ± 0.00
OPLS_2005	TIP4P-Ew	0.63	4.92 ± 0.17	0.50	1.27 ± 0.02	0.55	1.96 ± 0.04
AMBER ff99SB-ILDN	TIP3P	0.57	4.52 ± 0.09	0.39	0.78 ± 0.01	0.45	1.04 ± 0.01
AMBER ff03	TIP3P	0.49	1.53 ± 0.01	0.33	0.58 ± 0.00	0.36	0.65 ± 0.00
CHARMM27	TIP3P	0.53	9.03 ± 0.34	0.52	1.65 ± 0.03	0.44	1.00 ± 0.01
CHARMM22*	TIP3P	0.52	1.88 ± 0.02	0.44	1.03 ± 0.01	0.34	0.59 ± 0.00
OPLS_2005	TIP3P	0.59	11.65 ± 0.59	0.54	2.22 ± 0.04	0.55	3.37 ± 0.10
CHARMM27	mTIP3P	0.58	6.34 ± 0.20	0.47	1.20 ± 0.02	0.39	0.77 ± 0.01
CHARMM22*	mTIP3P	0.49	1.41 ± 0.01	0.50	0.76 ± 0.01	0.29	0.48 ± 0.00
OPLS_2005	TIP4P	0.60	8.18 ± 0.30			0.57	3.07 ± 0.07
AMBER ff03	TIP4P/2005	0.42	0.94 ± 0.01			0.29	0.46 ± 0.00
AMBER ff99SB-ILDN	SPC/E	0.60	2.98 ± 0.05			0.42	0.84 ± 0.01
CHARMM27	SPC/E	0.64	5.62 ± 0.17	0.49	1.23 ± 0.02	0.40	0.78 ± 0.01
CHARMM22*	SPC/E	0.51	1.54 ± 0.02			0.31	0.52 ± 0.00
OPLS_2005	SPC/E	0.66	5.26 ± 0.19			0.58	2.13 ± 0.04
experiment[22,45]c		∼0.25 ± 0.03	∼0.37 ± 0.05	∼0.22 ± 0.03	∼0.31 ± 0.05

Results are from 1-μs simulations and standard errors were calculated using a block averaging method.

Standard errors of Pbound were uniformly ≤0.01.

Experimental KA values of guanidinium and butylammonium acetate permit only a qualitative estimate of the associated error. Taking two experimentally measured KA values of guanidinium acetate using different protocols into account,[22,45] we estimate an error of ±0.05, although the true uncertainty is not known. Using this estimate, we have back-calculated the range of simulated Pbound values that would be expected in our simulation, based on the experimental KA.

Figure 2

Probabilities of binding (Pbound) between three different pairs of oppositely charged side-chain analogues using six biomolecular force fields with the TIP4P-Ew explicit water model. The Pbound values that correspond to the experimentally determined KA values of guanidinium acetate and butylammonium acetate are depicted as horizontal gray bars;[22] no experimentally measured KA is available for the imidazolium acetate system. Error bars represent 95% confidence intervals calculated using a block averaging method.[34]

One potential factor that could influence the degree of salt bridge formation in our simulations is the choice of force field parameters for the chloride ions. To verify that these parameters are not the cause for overestimating salt bridge strength, we carried out simulations of a single guanidinium/acetate pair (corresponding to concentrations of 0.1 M) with no chloride ions present. The resulting KA values are even higher than those measured in the presence of chloride ions, indicating that the chloride ion parameters do not cause disproportionate salt bridge stability (Table S3, Supporting Information).

Two obvious features of a force field that influence the strength of the salt bridges are the atomic charges and radii. As expected, the CHARMM22* force field yields KA values that are closer to experiment than those from the parent CHARMM27 force field since the atomic charges for the arginine, aspartate, and glutamate residues were parametrized specifically to reproduce the experimental association of guanidinium acetate.[3] However, the CHARMM22* force field does not produce as close agreement with experiment as the AMBER ff03 force field, which shows good agreement for butylammonium acetate. Given these results, it appears that the general strategy used to derive atomic charges for the AMBER ff03 force field is reasonably effective for modeling electrostatic interactions. This strategy involved the derivation of atomic charges in the presence of a continuum solvent model with a dielectric constant of 4 to mimic an organic solvent (protein-like) environment. The resulting atomic charges in the AMBER ff03 force field are notably less polarized than those in the AMBER ff9X family (including the AMBER ff99SB-ILDN force field tested here), which were derived in vacuum and share the same set of atomic radii. The AMBER ff03 atomic charges are also less polarized than those in the AMBER ff13α force field, with atomic charges possessing increased polarity relative to previous AMBER charge models.[4] In the AMBER ff13α charge model, nonpolarizable point charges have been fit to implicitly account for solvent polarization, using iterative cycles of classical MD simulations with explicit water (i.e., TIP4P-Ew) to estimate the water charge density around the solute, followed by quantum mechanical calculations to determine updated solute charges.

Interestingly, although certain critical atomic radii (e.g., the nitrogen in butylammonium and the oxygen in acetate) in the AMBER ff13α force field were adjusted from their original values in the AMBER ff99 force field to reproduce experimental hydration free energies of the relevant amino acid analogues,[4] the resulting strengths of the salt bridges are more overestimated, relative to the other tested AMBER force fields. Notably, the AMBER ff13α force field results in a free energy landscape for salt bridge formation that is significantly different from those of the other force fields. In particular, as shown by the PMFs as a function of the minimum N–O distance between the oppositely charged analogues for the three types of salt bridges (Figure 1), the free energy minima for the bound states are consistently shifted to the right in the AMBER ff13α force field relative to the other force fields. When we substituted the atomic radii in the AMBER ff13α force field with the original radii from the AMBER ff99 force field, the free energy minima for the bound states shifted back toward those of the other force fields and yielded significantly deeper minima as well as more pronounced desolvation barriers, particularly for the guanidinium/acetate and imidazolium/acetate systems (Figure S2, Supporting Information).

Since the atomic charges of the OPLS_2005 force field are not significantly different from those of the other force fields, the most likely reason for the fact that this force field overestimates the KA values to the greatest extent is that the atomic radii of the nitrogen-attached hydrogen atoms are smaller than those used by the other force fields, potentially allowing the pairs to associate more closely and increasing their electrostatic attraction (Table S1, Supporting Information). Consistent with this notion, simulations using the OPLS-AA force field, which differs from the OPLS_2005 force field only in that it omits atomic radii for these hydrogen atoms, resulted in slightly more strongly associated salt bridges (Table S2, Supporting Information). We note that the ranking of the strengths of the three types of side-chain salt bridges in our study by the AMBER ff99SB-ILDN, CHARMM27, and OPLS-AA force fields is consistent with that observed for their oppositely charged termini in a recent study by others.[21]

Association Constants of Amino Acid Dipeptides

As mentioned above, we additionally tested salt bridge formation of the Arg/Asp pair by simulating association/dissociation of blocked amino acid dipeptides, testing six different force fields in conjunction with the TIP4P-Ew explicit water model. As shown in Table 2, the relative ranking of the force fields in terms of the KA is generally consistent with our results from the corresponding side-chain analogue system (guanidinium/acetate). The only exception is the AMBER ff13α force field, which yields the weakest KA for the association of the amino acid dipeptides, as opposed to an intermediate KA value for the association of guanidinium/acetate. As indicated by the PMF between the arginine and aspartate dipeptides (Figure 3), the bound state free energy minimum of the AMBER ff13α force field is the most shallow among the tested force fields, corresponding to the lowest frequency of salt bridge formation. The inclusion of the backbone groups, therefore, appears to alter its propensity for salt bridge formation, likely through the competition of side-chain/backbone interactions with the side-chain/side-chain interactions between the two amino acids. This result emphasizes the benefit of using unbiased simulations; had the relative orientations of the amino acids been fixed as in previous studies,[14] any effects of significant side-chain/backbone interactions on the frequency of salt bridge formation would not have been apparent.

Table 2

Association Constants (KA), Probabilities of Side-Chain/Side-Chain (PboundSC/SC), and Side-Chain/Backbone Association (PboundSC/BB), Respectively, for Blocked Arginine and Aspartate Dipeptides Using Six Biomolecular Force Fields with the TIP4P-Ew Water Modela

		arginine/aspartate
force field	water model	KA (M–1)	PboundSC/SC	PboundSC/BB	PboundSC/SC/PboundSC/BB
AMBER ff99SB-ILDN	TIP4P-Ew	3.5 ± 0.3	0.041 ± 0.003	0.010 ± 0.001	4.1 ± 0.4
AMBER ff03	TIP4P-Ew	2.8 ± 0.2	0.033 ± 0.002	0.017 ± 0.001	1.9 ± 0.2
AMBER ff13α	TIP4P-Ew	1.9 ± 0.1	0.023 ± 0.002	0.018 ± 0.001	1.3 ± 0.1
CHARMM27/CMAP	TIP4P-Ew	11.4 ± 1.2	0.113 ± 0.012	0.007 ± 0.001	16.4 ± 2.0
CHARMM22*	TIP4P-Ew	3.3 ± 0.2	0.038 ± 0.003	0.007 ± 0.001	5.5 ± 0.5
OPLS_2005	TIP4P-Ew	17.0 ± 1.9	0.153 ± 0.016	0.013 ± 0.001	11.8 ± 1.6

Results are from 10-μs simulations and standard errors were calculated using a block averaging method.

Figure 3

Potentials of mean force (PMF) between blocked arginine and aspartate dipeptides using six biomolecular force fields with the TIP4P-Ew water model. The larger noise level compared to the data presented in Figure 1 is caused by simulating a single pair of binding partners, rather than a concentrated solution.

To monitor side-chain/backbone association, we used the same minimum N–O distance coordinate and bound and unbound state definitions as used for the side-chain/side-chain interactions. A comparison of the relative probabilities of side-chain/side-chain versus side-chain/backbone association (Figure 4) reveals that the force fields generally prefer side-chain/side-chain association by a factor of ∼2 or more over side-chain/backbone association. The exception is AMBER ff13α, which shows a lower preference for side-chain/side-chain association of ∼1.3. This slight preference over side-chain/backbone association is likely due to the substantially more polarized backbone amide and carbonyl groups of the AMBER ff13α force field relative to previous AMBER force fields (including the AMBER ff99SB-ILDN and AMBER ff03 force fields).[4] Thus, as a result of the delicate balance of side-chain/side-chain and side-chain/backbone interactions, the strength of the Arg/Asp salt bridge appears to be most accurately modeled (least overstabilized) by the AMBER ff13α force field in a model system that is representative of a protein environment.

Figure 4

The probability of side-chain/side-chain association (PboundSC/SC) over the probability of side-chain/backbone association (PboundSC/BB) for blocked arginine and aspartate dipeptides using six biomolecular force fields with the TIP4P-Ew water model. Error bars represent 95% confidence intervals calculated using a block averaging method.[34]

Influence of the Water Model

In addition to the force field, the choice of water model can affect the strength of salt bridges. To evaluate the influence of the water model, we tested the above three side-chain analogue systems with a selection of force field/water model combinations in addition to the force field/TIP4P-Ew combinations. Regardless of the water model, the relative ranking of the force fields is unchanged in terms of the KA values, with Pbound varying by 5–10% between the water models (Table 1). We also evaluated the dependence of salt bridge interactions on the dielectric constant of the employed water model (εwater). Interestingly, despite the fact that the SPC/E water model yields a computed dielectric constant (εwater = 70; Table S4, Supporting Information) that is closest to the experimental value (εwater = 78.4)[46] among all of the tested water models, the use of the SPC/E water model results in stronger salt bridge interactions than seen with the TIP4P-Ew water model, in which the εwater value is underestimated (εwater = 56; Table S4, Supporting Information). In fact, as shown in Figure 5, there appears to be no clear correlation between εwater of the water model and the strength of the salt bridges. As an aside, the CHARMM27 and CHARMM22* force fields were tested with both the standard TIP3P and the CHARMM-modified TIP3P (mTIP3P) water models with which they were developed. The mTIP3P water model includes atomic radii on hydrogen as well as oxygen atoms,[42] whereas standard TIP3P includes only the oxygen atom. Using the mTIP3P water model consistently results in a lower KA values, in better agreement with experiment. This suggests that it may be advisible to use of the CHARMM-modified TIP3P, rather than the standard TIP3P water model with any CHARMM force field.

Figure 5

Relationship between the probabilities of binding (Pbound) and the dielectric constant of the water model (εwater) for three different pairs of oppositely charged side-chain analogues. The displayed data are simulation results using six biomolecular force fields and six explicit water models. The Pbound values that correspond to the experimental association constants (KA) for the guanidinium acetate and butylammonium acetate systems are depicted as horizontal gray bars;[22] no experimentally measured KA is available for the imidazolium acetate system. The εwater values were calculated from the first 100 ns of each simulation. For each model, the presence of the solutes lowers the εwater for each system by 10–15 relative to that of pure water (Table S4, Supporting Information). Note that the error bars are not visible since 95% confidence intervals for both Pbound and εwater lie within in the symbols’ area in the graph.

Conclusions

We compared the modeling of salt bridge interactions using six current biomolecular force fields. Three different salt bridges (Arg/Asp, Lys/Asp, and His(+)/Asp) were simulated and considerable differences in their strengths were noted, both between the force fields and with experiment. Given the availability of experimentally measured KA values for the association of oppositely charged side-chain analogues, we have focused primarily on modeling salt bridge formation using these systems. We also tested the applicability of our results to amino acids by simulating blocked amino acid dipeptides for one of the salt bridges, Arg/Asp.

Our side-chain analogue simulations reveal that the computed KA values are generally overestimated, relative to experimental values, with the AMBER ff03 force field overestimating the strengths of the salt bridges to the least extent and the OPLS_2005 force field to the greatest extent when using the same water model (TIP4P-Ew). For the blocked arginine and aspartate dipeptides, we observed general agreement in the relative ranking of the force fields with that obtained from simulations with the corresponding side-chain analogues. The only exception is the AMBER ff13α force field, which resulted in the lowest probability of salt bridge formation, likely due to the presence of competing side-chain/backbone interactions. Thus, while the AMBER ff03 force field overestimates the strength of salt bridges to the least extent for the side-chain analogue systems, the AMBER ff13α force field results in an even lower frequency of salt bridge formation than the AMBER ff03 force field for the complete amino acids. Finally, we examined the influence of the water model on the strengths of the salt bridges. Irrespective of the water model, the relative ranking of the force fields remained unchanged, with no clear correlation between the probability of binding (salt bridge formation) and the dielectric constant of the solvent (εwater).

In conclusion, when running MD simulations in which salt bridge formation may be of interest, careful attention should be paid to the specific force field and water model in simulations of protein systems. Several current force fields yield considerably higher KA values than those experimentally determined, a discrepancy that may lead to erroneous conclusions. Our encouraging results with the AMBER ff13α force field suggest that charge derivation strategies that implicitly incorporate solvent polarization from explicit water may significantly extend the lifetime of fixed-charge force fields, which include all of the force fields tested in this study. Nonetheless, departures from these simple point charge models may also be necessary. For example, using polarizable force fields that permit varying the charge distribution within a molecule based on both its conformation and environment may alleviate such shortcomings. In the past, the solvation of ions[47,48] and charged small molecules,[49] have been modeled using polarizable force fields, resulting in improved agreement with experiment, compared to CHARMM27 and AMBER ff99 force fields (equivalent to the AMBER ff99SB-ILDN force field used here). Future work will determine whether or not this also holds for protein salt bridges.