How well do force fields capture the strength of salt bridges in proteins?

Introduction

Proteins are stabilized via the concerted action of numerous weak forces including those that arise from hydrogen bonds, the hydrophobic effect and salt bridges (Dill, 1990; Zhou & Pang, 2017). While we now know much about the relative contributions and physical origins of these effects, we still do not have quantitative models that allow us, for example, to predict accurately the overall stability of a protein given its three-dimensional structure. A quantitative understanding of protein stability would aid both to design proteins with improved stabilities (Foit et al., 2009), as well as understand how loss of protein stability may give rise to disease (Casadio et al., 2011; Nielsen et al., 2017).

We here focus our attention on salt bridges in proteins, where they occur between oppositely charged ionisable groups, most commonly between the negatively charged side chains of aspartate and glutamate or the carboxy-terminus, and positive charges in arginine, lysine or histidine side chains, or with the N-terminal ammonium group. A salt bridge is formed when the two residues involved are spatially close enough to form energetically favourable electrostatic interaction between the (partially) charged atoms (Barlow & Thornton, 1983). Although the first quantitative models of protein electrostatics are almost one hundred years old (Linderstrøm-Lang, 1924), our ability to predict electrostatic properties remain incomplete. Indeed, while ionic interactions were also early suggested to contribute substantially to the stability of certain proteins (Speakman & Hirst, 1931), their importance and energetic contribution was already controversial before they had been observed in experimentally-derived protein structures (Jacobsen & Linderstrøm-Lang, 1949).

The stability of a salt bridge depends on the environment around the residues, the pH, and the distance and geometric orientation between the involved residues. Salt bridges may be located on the protein surface, where they are exposed to the solvent, or they can be found buried within the hydrophobic interior of the folded protein. Studies on different protein families have found that buried salt-bridges are more likely to be conserved and functionally relevant than surface exposed ones (Schueler & Margalit, 1995; Takano et al., 2000). Surface exposed salt bridges are generally weaker, more variable in their contribution to stability, and more difficult to predict (Sarakatsannis & Duan, 2005). Double mutant cycles in barnase illustrate the distinction between forming a salt bridge in the low dielectric environment of the protein core versus the higher dielectric on the protein surface. These experiments show that an internal salt bridge stabilizes the protein by more than 3 kcal mol⁻¹ (Vaughan et al., 2002), whereas a surface exposed salt bridge is about 10 times weaker, ∼0.3 kcal mol⁻¹ at low ionic strength, a value that decreases to zero at higher ionic strength (Serrano et al., 1990).

Although solvent-exposed salt bridges are often observed in crystal structures of proteins, it is thus not always clear whether these interactions are stable in solution, in particular because their dynamic and transient nature make them difficult to study. With its site-specific resolution and ability to detect even transient interactions, NMR spectroscopy can, however, be used to study salt bridges in solution, and to determine the extent to which different residues interact. In one intriguing study, a range of different NMR experiments were used to examine three different potential salt bridges in the B1 domain of Protein G (hereafter called GB1) (Tomlinson et al., 2009). Six different crystal structures of GB1 show that three of six solvent-exposed lysines, K12, K39 and K58, form salt-bridges to the nearby acidic residues E12, E23 and D55 (Tomlinson et al., 2009) (Fig. 1A). To examine whether these ionic interactions are present in solution, the authors performed a series of experiments that include monitoring both the lysine nitrogen and proton chemical shifts, and the hydrogen-deuterium isotope effects on the ammonium group, while titrating the carboxylates to protonate them. Surprisingly, for two of the putative ion pairs (K12–E23 and K58–D55) the results showed essentially no observable changes at the lysines as the carboxylates changed protonation state, suggesting no substantial formation of a salt bridge. For one pair, the intra-helical K39–E35, there was a small change in the ¹⁵N chemical shift of the ammonium group as pH was varied to protonate E35 (and other carboxylates), though no change was observed for the ¹H chemical shift nor for the isotope effects. These observations were also supported by the determination of the pKa values of the carboxylate and ammonium groups in these residues. The authors therefore conclude that two of the salt bridges (K12–E23 and K58–D55) are not formed in solution, while the third (K39–E35) may be weakly formed (though likely at a low population).

Computational methods provide an alternative approach to study the structure and dynamics of ionic interactions in solution (Kumar & Nussinov, 2002). All-atom, explicit solvent molecular dynamics (MD) simulations, in particular, can be used to provide insight into the structure and energetics of biomolecules in an aqueous environment. For such simulations to provide an accurate description of salt bridges it is, however, important that the energy function (force field) used provides an accurate description of the balance between the many forces that determine salt bridge strengths. As part of the re-parameterization and validation of one of the more accurate force fields, CHARMM22*, we discovered that ionic interaction between guanidinium and acetate ions were ∼10 times too strong in the original CHARMM22 force field (Piana, Lindorff-Larsen & Shaw, 2011). This observation led to a re-parameterization (the “DER correction”) of the partial charges in Asp, Glu and Arg side chains, to bring the interactions closer to experiment (but still slightly overestimated) (Piana, Lindorff-Larsen & Shaw, 2011). Building upon these ideas, Debiec and colleagues examined the interactions between analogues of both arginine, lysine and histidine side chains with carboxylates in a range of force field and water models (Debiec, Gronenborn & Chong, 2014). Although they found considerable variability between different force fields, the general observation from the CHARMM force fields hold true, namely that the ionic interactions between side chain analogues are overstabilized relative to experiment. Such overstabilization is not without consequence when studying protein dynamics using simulations. For example, in a study of voltage gating in potassium channels the simulations were performed with the DER correction, or the even more substantial DER2 correction, in order to be able to observe the switch in ionic interactions that are central in channel gating (Jensen et al., 2012). Similarly, when using free energy perturbation MD simulations to predict the change in protein stability from a mutation, the results are less accurate when the mutations involve a change in charge (Steinbrecher et al., 2017).

In addition to the DER correction described above, several groups have examined and modified the charges in protein force fields in order to provide a more accurate description of protein electrostatics. For example, Jensen (2008) updated the parameters for the charged amino acids in the OPLS force field to obtain a self-consistent set of parameters that reproduce experimental hydration free energies for side chain analogues. Debiec et al. (2016) used the Implicitly Polarized Charge (IPolQ) method to derive a new set of partial charges for Amber, resulting in the ff15ipq force field, and showed that these modifications substantially improved the agreement with experiments on the association between analogues of charged side chains.

Here we build on the work described above by examining the extent to which seven different fixed charge force fields capture the behaviour of the three putative salt bridges in GB1. In particular, we conducted explicit-solvent molecular simulations of GB1 using the Amber ff99SB*-ILDN, Amber ff03w, CHARMM27, CHARMM22*, OPLS_2005, AMBER ff15ipq and a99SB-disp force fields. In contrast to studies of side chain analogues, simulations of the salt bridges in the context of a folded protein takes into account the natural geometrical and energetic constraints imposed by the protein scaffold. To accelerate sampling and ensure convergence we used metadynamics simulations to map the free energy surface of salt bridge formation. The results reveal that also in the context of a folded protein most of these force fields overestimate the population of ionic interactions, and also provide insight into the geometry of salt bridges in proteins.

Materials and Methods

In order to investigate potential salt-bridging interactions between residues K12-E23, K39-E35 and K58-D55 in GB1 we performed both standard MD simulations as well as enhanced sampling simulations using well-tempered metadynamics (WT-MetaD) (Barducci, Bussi & Parrinello, 2008). To test the effect of the force field we used seven different force fields: Amber ff99SB*-ILDN (Hornak et al., 2006; Best & Hummer, 2009; Lindorff-Larsen et al., 2010), Amber ff03w (Best & Mittal, 2010), CHARMM27 (CHARMM22 with the CMAP correction) (MacKerell et al., 1998; Mackerell, Feig & Brooks, 2004), CHARMM22* (Piana, Lindorff-Larsen & Shaw, 2011), OPLS_2005 (Banks et al., 2005), a99SB-disp (Robustelli, Piana & Shaw, 2018) and AMBER ff15ipq (Debiec et al., 2016) (Table 1). Each force field was combined with its “native” water model, namely the CHARMM-specific TIP3P water model for CHARMM force fields (MacKerell et al., 1998) (herein called TIPS3P), the original TIP3P water model (Jorgensen et al., 1983) for the Amber ff99SB*-ILDN force fields, TIP4P/2005 (Abascal & Vega, 2005) for Amber ff03w and the SPC/E water model (Berendsen, Grigera & Straatsma, 1987) for OPLS. In the case of the AMBER ff15ipq force field, we used the SPC/E_b water model, and the simulation box was constructed with AmberTools 17 (Case et al., 2017) and then converted to Gromacs format. A slightly modified version of the TIP4P-D water model (Piana et al., 2015) was used for the a99SB-disp force field (Robustelli, Piana & Shaw, 2018). We also examined the effect of the water model by testing CHARMM22* with standard TIP3P.

Results

Molecular dynamics simulations of salt bridge formation

We performed MD simulations of GB1 using eight different combinations of force fields and water models (Table 1) to examine the formation of salt bridges, and to benchmark against NMR experiments. While unbiased MD simulations showed reversible formation and breaking of the ionic interactions (Fig. 2A), we decided to use enhanced sampling metadynamics simulations (Fig. 2B) to ensure better convergence of the free energy of salt bridge formation. In both types of simulations, we find that the carboxylate and ammonium group spend time at multiple, relatively distinct sets of distances, though with differences between force fields and salt bridges (see below). Thus, for each of the eight different force field combinations we performed three metadynamics simulations, one for each salt bridge pair. In each of these 24 simulations, we simultaneously biased and monitored two CVs, corresponding to the distance between the ammonium nitrogen atom in the lysine side chain, and each of the two oxygen atoms in the carboxylate group in aspartate or glutamate. The metadynamics bias did not have any major effect on the protein stability over the course of the 100 ns simulations. For example, the average backbone RMSD was below 1.5 Å for both the unbiased and biased simulations (Fig. S1).

Assessing convergence of simulations and effect of simulation parameters

When MD simulations are in apparent disagreement to experiments, it is important to assess whether these differences are due to the force field and simulation parameters, or whether they can be explained e.g., by insufficient sampling. While it is difficult to prove that simulations are converged, we performed several tests that suggest that our observations are rather robust. First, as discussed above, the high degree of symmetry of the calculated free energy surfaces suggest that the individual nitrogen–oxygen distance distributions are relatively converged. Second, we calculated the free energy of salt bridge formation during the simulation (from the populations where the nitrogen–oxygen distance is <5.0 Å) and monitored the time evolution of this (Fig. 3 shows the results from three WT-MetaD simulations with CHARMM22*/TIPS3P). While the values fluctuate for the different salt bridges, they are relatively stable in the second half of the simulations. Finally, for the K35–E39 salt bridge we also performed a parallel-tempering metadynamics simulation in CHARMM22*/TIPS3P, and compared the results to that from the WT-MetaD (Fig. 3). The two simulations converge to essentially the same free energy difference. Based on these observations we estimate an error of about 0.5 kcal mol⁻¹, corresponding roughly to an error of 0.2 for the population of the salt bridges when the free energy of salt bridge formation is close to zero.

Another key parameter in a simulation is the cutoff used to truncate the Lennard-Jones interactions and to switch between a direct calculation of electrostatic forces and the calculations of longer-range electrostatics using Ewald summation. For the K12–E23 and K35–E39 salt bridges and the CHARMM22*/TIPS3P force field we therefore repeated the WT-MetaD simulations varying this cutoff between 9.0 Å–14.0 Å and calculated the free energy of salt bridge formation (Fig. 4). While the results vary slightly and in accordance with the estimated uncertainty of our simulations, we find no systematic dependency of the free energy differences on the cutoff used. Thus, these calculations suggest that the cutoff used in the simulations are sufficient to obtain reasonable results.

Discussion

We have used enhance sampling simulations to analyse the stability of salt bridges in commonly used fixed charge force fields for molecular dynamics simulations. We estimate the error of the calculated free energies of salt bridge formation to be ∼0.5 kcal mol⁻¹, corresponding to errors of the populations of ∼20% (down to 10% for the most skewed populations). Surprisingly, we find that the K12–E23 and K58–D55 pairs form the most stable salt bridges in most of the force fields (∼64% on average across force fields) with the K39–E35 pair being distinctly less stable (on average ∼35%). These observations are in clear deviation from those expected from the experimental observations, in that the simulations generally appear both to overestimate the stability of the salt bridges, and that the order of salt bridge stability also appears to be wrongly predicted. In addition to differences in the partial charges we note that differences in bonded terms, e.g., in the torsional terms that affect side chain rotamers, can affect the relative stability and geometry of the salt bridges.

The Amber ff15ipq force field is the major outlier to the above observed trends. Thus, in contrast to the other force fields that overestimate the stability of the salt bridges, this force field is in substantially better agreement with experiment in that it finds that the K12–E23 and K58–D55 salt bridges are very weakly formed (10% and 2%, respectively). The K39–E35 salt bridge is formed ∼20% of the time in the ff15ipq force field, and a more detailed comparison to the raw NMR data would be needed to determine whether this is in agreement with experiments.

For the intra-helical K39–E35 interaction we note that K39 also forms transient interactions with other residues (including E48 and N43), suggesting that the variations observed between force fields might be the cumulative effect of a number of differences between the force fields. We find that the ionic interactions are slightly stronger in TIP3PS than in standard TIP3P, in line with previous observations from model compounds (Debiec, Gronenborn & Chong, 2014). As also observed before (Piana, Lindorff-Larsen & Shaw, 2011; Debiec, Gronenborn & Chong, 2014), we find that the DER correction in CHARMM22* decreases the strength of salt bridges to be more in line with the experiment.

Conclusions

Ionic salt bridge interactions are pervasive in experimental protein structures and have in certain cases been shown to contribute both to protein stability (Vaughan et al., 2002) and fast association kinetics (Schreiber, Haran & Zhou, 2009). Nevertheless, our understanding of the geometry and energetics of salt bridges in solution is limited by the difficulty in experimental and computational studies. While molecular dynamics simulations in explicit solvent may in principle be a quantitative and predictive model for salt bridge formation, the accuracy of such simulations hinges upon the force fields used. Building upon earlier work on model compounds (Piana, Lindorff-Larsen & Shaw, 2011; Debiec, Gronenborn & Chong, 2014; Debiec et al., 2016), we have here performed simulations of GB1 to determine the free energy landscape of lysine-carboxylate salt bridge formation. In these analyses we opted to examine only fixed charge force fields, but note that at least in the case of guanidinium-acetate interactions polarizable force fields have been found to give good agreement with experimental data on salt bridge strengths in model compounds (Debiec et al., 2016). Our comparison with experimental NMR data on GB1 suggest that while the force fields recapitulate the transient and weak nature of these solvent exposed ionic interactions, all but one appear to overestimate their stability slightly. This observation is in line with those from the model compounds, suggesting that together these kinds of calculations might also be used to improve the force fields. Indeed, the newly developed charge model for the Amber ff15ipq force field provides a more balanced description of salt bridge interactions for model compounds (Debiec et al., 2016), and our results on GB1 support this observation.

In general, we urge practitioners of MD simulations to take the small, but significant force field bias in many force fields into account when interpreting the importance of salt bridges observed in simulations, unless such observations are supported by experimental data. Finally, we hope that experimentalists will continue to develop approaches to study electrostatic interactions in proteins, in particular experiments that can be compared directly to simulations. Recent examples include extension of the GB1 studies to salt bridges in barnase (Williamson et al., 2013), novel NMR methods for studying electrostatics (Hass & Mulder, 2015), NMR methods to study arginine side chains (Mackenzie & Hansen, 2017; Yoshimura et al., 2017), and the engineering of a protein without titratable side chains as a platform for studies of protein electrostatics (Højgaard et al., 2016).

Supplemental Information