A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization

Calcium signaling affects a plethora of cellular processes with specific proteins designed by nature as Ca²⁺ binders.¹ A necessary condition for calcium to act as a signaling agent¹ is that the interaction with its binding partner has to be weak enough (i.e., below ten k_BT) so that it can bind and unbind at biologically relevant timescales.² This is an imperative condition that needs to be satisfied by any molecular simulation that aspires to faithfully model processes involving calcium signaling. However, the widely used classical non-polarizable force fields tend to dramatically overbind calcium and other ions of high charge density (such as magnesium or lithium) not only to proteins and peptides^3–6 but also to phospholipids^7–11 and nucleic acids.^12,13 Similar overbinding is observed for other ions too,¹⁴ and the issues are further manifested in the overestimated salt bridge strength between charged protein residues.^15,16 Attempts to cure the overbinding problem by introducing an ad hoc steric repulsion term between the binding partners may improve the situation but provide, at best, the right result for an often unphysical reason. Physically, the overbinding originates from the lack of electronic polarization. The fast electrons should always be at equilibrium with nuclei positions and, hence, effectively screen the ion–ion Coulomb interaction.¹⁷ This missing screening causes the observed disbalance between ion–water and ion–counterion interactions in non-polarizable force fields, which can, in principle, be fixed by explicitly introducing electronic polarizability. This way, one might reproduce the realistic behavior by incorporating missing physics in the model. However, due to issues connected with parameterization, numerical stability (i.e., the so called polarization catastrophe¹⁸), and increased computational costs, biomolecular simulations with explicitly polarizable force fields are still rare.^19–21 

Here, we address the question whether one can get a satisfactorily accurate behavior for (almost) correct reasons without extra computational cost. As demonstrated below, the answer is yes, namely, by including electronic polarization effects in a mean-field way via charge scaling. In this Perspective, we aim to give an overview of recent efforts and strategies toward the development of scaled-charge force fields within the electronic continuum correction (ECC) approach, with a focus on applications in systems of biological interest. We want to provide a cookbook on how one can develop, test, and use such force fields as well as highlight the potential pitfalls, open questions, and limitations of such an approach. We will first review the various existing strategies used to correct ion binding artifacts and ion pairing in simulations, including the ECC strategy. We will then present the different steps toward the application of the ECC framework to increasingly complex biosystems, starting from simple salt solutions all the way up to protein- and lipid-containing systems. While we will do our best to mention all contributions to this field, a particular focus will be put on the work performed in our group.

One strategy to improve the aforementioned overbinding artifacts is to explicitly include electronic polarization into classical force fields. This can be done in different ways, the most common approaches being the induced point dipole, the fluctuating charge, and the Drude oscillator models. These approaches and their applications in the simulation of electrolytes, ionic liquids, or biomolecular systems have been reviewed elsewhere,^22–24 and a detailed description thereof is outside the scope of the present review. Polarizable force fields have allowed for more realistic simulations of ionic liquids and electrolytes.^24,25 Focusing on the interactions of ions with biomolecules, they have been successfully used to investigate, for instance, the geometry and relative stability of Mg²⁺ binding modes to phosphate-based species^26–28 and the cation specificity in proteins.²⁹ Polarizable force fields have also been shown to better capture the strength of salt bridges.¹⁶ 

However, despite these significant advances and recent developments, polarizable force fields are not yet broadly used by the biosimulation community. The vast majority of simulations of complex biomolecular systems are still performed with non-polarizable force fields, even if they suffer from well-identified artifacts. The reasons are that due to a higher number of parameters, complex molecules are, in general, tedious to parameterize for polarizable force fields, which often require the use of specialized software with a smaller user base. Additionally, their slower computational performance hampers their application to large biomolecular systems. Fortunately, significant advances have recently been made to overcome some of these limitations.^30,31

The current limitations in the explicit treatment of electronic polarizability forced the simulation community to look for simpler, easy-to-use, and readily available methods to improve ion binding and pairing in widely used simulations. Several approaches have thus been developed to compensate for the artifacts arising from the lack of electronic polarization in such simulations within the existing interaction terms of non-polarizable force fields. In the following, we will briefly mention the main strategies adopted in the literature, before focusing more closely on the electronic continuum correction approach.

One strategy to take into account induced polarization in a mean-field approach has been to bypass the commonly used combination rules for Lennard-Jones parameters. Ad hoc pair-specific Lennard-Jones parameters (for example, the NBFIX corrections in the CHARMM force field³²) are defined for specific interacting pairs of ions or charged groups in biosystems^33–36 in order to improve their pairing properties.^7,36–38 There are several potential limitations to these approaches. The parameterization of each charge pair of interest has to be done independently, which can be tedious. Furthermore, experimental data required as parameterization targets are not always available and are often difficult to obtain. This is, of course, a problem for any parameterization effort, but it can often be partially overcome by relying on combination rules instead of parameterizing each interaction pair separately. The result is that the efforts relying on pair-specific interaction terms only focus on the interactions between a few atom types at a time. It is also unclear how the introduction of these specific parameters affects the balance of the interactions of other (unmodified) pairs and whether the corrections developed in different studies are mutually compatible and transferable to other force fields. Finally, in situations where the binding of a charged moiety is critical, the strategy of simply pushing away charges to decrease their interaction results in an inaccurate description of the interactions. Both polycoordination, e.g., Ca²⁺ binding sites in proteins, and interactions where directionality is fundamental, e.g., Zn²⁺ binding sites, are problematic.

Other authors suggested³⁹ adding a 1/r⁴ term to the standard 12-6 Lennard-Jones potential to account for charge-induced dipole interactions that are neglected in standard force fields.³⁹ A force field with the resulting 12-6-4 Lennard-Jones-type non-bonded model has been developed for a series of monovalent⁴⁰ and multivalent⁴¹ ions in combination with different water models. This model was shown to be able to correctly capture hydration free energies, ion–water oxygen distances, and coordination numbers simultaneously and has given promising results on ion pairing properties. However, for divalent cations such as Mg²⁺ or Zn²⁺, further tuning of their pairwise interaction with specific RNA groups, e.g., its phosphate backbone or the guanosine and adenosine N7 atom, was found necessary to reproduce experimental binding free energies.⁴² 

Instead of modifying the parameters of the van der Waals dispersion interaction, different approaches have suggested to directly tune the electrostatic properties of the solute through the modification of atomic charges. In this spirit, refinements^43,44 of standard force fields have been suggested to reduce the strength of protein salt bridges that are strongly overestimated in most standard force fields.^15,16 In particular, the ff15ipq force field,⁴⁵ recently developed with the implicitly polarized charge method (IPolQ), uses non-polarizable atomic charges that were rederived to implicitly capture the energy of polarization due to the presence of a solvent, namely, water in biological systems. This force field, which keeps full integer charges on charged residues but changes the charge distribution compared to the original force field, has been shown to significantly improve the description of the strength of salt bridges.

Another approach is to scale the total charge of ions or charged moieties to improve their pairing properties. Charge scaling has been for a long time prevalent in the context of ionic liquid simulations,^24,46,47 where a scaling factor of 0.7–0.8 has been shown empirically to give reasonable thermodynamic properties.⁴⁸ Charge scaling is not an alien concept in biological simulations either,^49,50 though its adaptation has been so far scarce and mostly limited to issues involving ion channels⁵¹ or enzymatic catalysis.⁵² In recent years, the charge scaling strategy has been revived and physically grounded in a series of papers by Leontyev and Stuchebrukhov^53–56 who rationalized this approach in terms of a mean-field theory called here the electronic continuum correction, which we will focus on from here on.

A physically well-justified approach to decrease the binding between charged moieties and groups in non-polarizable force field simulations is to introduce electronic polarizability within a mean-field approach via ECC. Note that the missing electronic polarizability in non-polarizable force fields leads to insufficient screening of interactions between charges. As the electronic polarizability is fairly constant in biological systems, one can incorporate it as a dielectric continuum, which can, in principle, be done in two ways—either by immersing the system in a continuum medium with the electronic polarizability equal to ϵ_el (ECCϵ)⁵⁷ or by scaling charges of an uncorrected force field. Both approaches are theoretically equivalent. In practice, the uniform weakening of all Coulombic interaction by a background ϵ_el is unpractical for biomolecular simulations, as such a global scaling of the charges may introduce undesired changes (overscaling) in carefully parameterized interactions. Features such as dihedral angular terms in molecules, sensitive to changes in 1–4 charge interactions, or the area per lipid in a membrane, sensitive to modifications in the overall lipid–lipid interactions, are often carefully fitted in standard force fields to reproduce experimental data and should thus be left untouched. Alternatively, one may develop an entirely new force field, but this is an enormous task considering the complexity of biomolecules and their environments. ECCϵ has been applied with promising results for some simple systems⁵⁷ and seems to resolve some issues related to charge scaling.⁵⁸ However, the charge scaling ECC variant is less intrusive, as it relies on the fact that only heavily charged groups or atoms contribute significantly to charge–charge interactions. This means that partial charges on weakly polarized groups can be safely left untouched.

The possibility to effectively account for electronic polarization via charge scaling has been around since the 1980s. A pioneering work by Jönsson, Edholm, and Teleman stated explicitly: “One very simple way to remedy the lack of dielectric screening in the water potential is to scale ionic interactions with the high frequency dielectric constant….”⁴⁹ However, later, in the same paper, the authors scaled the ionic charges by a factor of 0.5, which from the present point of view seems like overscaling. This ambiguity may explain why systematic charge scaling has not caught up until about a decade ago. Only then the appropriate scaling factor for ionic charges—equal to the inverse square root of the high-frequency dielectric constant of water (i.e., 0.75)—was introduced.⁵³ In the meantime, nevertheless, charge scaling has been applied as an empirical correction to several problems, including studies on the structural and dynamic properties of ionic liquids,⁵⁹ charge transfer from catalytic ions in enzymes,⁶⁰ and ion channels.⁵¹ 

The seminal contribution of Leontyev and Stuchebrukhov lies in developing charge scaling on a firm physical basis. By a simple manipulation of Coulomb’s law, one can show that scaling the charges by the inverse square root of the high frequency dielectric constant is mathematically equivalent to immersing the charges into a medium with the respective relative permittivity.^54–56,61 Indeed, following the derivation of classical electrostatics, one can show that such charge scaling is rigorous and represents a way to preserve the symmetry between the “source” and “test” charges in a dielectric medium.⁶² This mean-field approach is physically well justified for systems that are roughly homogeneous in their electronic responses, i.e., which do not possess spatial discontinuities in their high-frequency dielectric constants. Fortunately, most of the biologically relevant molecular environments composed of moieties such as water, proteins, lipids, and nucleic acids are quasi-homogeneous in this respect with a high-frequency dielectric constant of about 2 within the relevant temperature range (see Table I), which makes the ECC charge scaling strategy applicable to a wide range of systems, including some of those containing interfaces.

Strictly speaking, the applicability of the ECC approach is limited to systems with an essentially uniform high-frequency dielectric constant across the system.⁶³ Therefore, systems where an abrupt change of the high-frequency dielectric constant occurs and where scaled charges can enter such interfacial regions are strictly speaking not compatible with the ECC approach. Fortunately, in most of the biologically relevant interfaces, at least one of these two conditions is not met. Still, simulations of a solvent–air or solvent–metal interface are problematic.⁶⁴ In such cases, on the solvent side, the ECC scaling factor should be applied, whereas on the other side, a different scaling factor should be applied, e.g., a scaling factor of 1 (full charges) for air. Interestingly though, the application of ECC in molecular dynamics simulations of solid–liquid interfaces disregarding the limitation seems to improve the interaction of charges with highly polar surfaces.^65–69 Such changes are not currently handled by any MD simulation codes.

In principle, the electronic continuum correction (ECC) accounts for electronic polarization in a mean-field way as follows: (i) Scale any charged group so that its total charge is 0.75 times its nominal charge. This takes care of electronic polarization by effectively putting the charged moiety in an environment with a high-frequency dielectric constant of water ϵ_el ∼ 1.78 (corresponding to a refraction index of $n=ϵel∼1.34$⁠). The charge scaling factor is then $1ϵel=0.75$⁠, which is to a good approximation valid also for any aqueous environment or complex biological system.⁷⁰ (ii) Run a standard classical molecular dynamics that recovers naturally the nuclear polarization.^53–56 At least, this is how it should be done in the ideal world where the water model itself is fully compatible with the ECC approach, i.e., its dielectric constant only contains the nuclear polarization with ϵ_nuc = 45. The existence of such an ideally behaving ECC water model is not guaranteed, but a desirable target for the ECC ecosystem. Until such a water model becomes available, the ECC development will rely on readily available non-ECC water models that have provided undoubted success in distinct molecular environments.^{4–6,10,63,71,72} The commonly used non-polarizable water models exhibit dielectric constants between 50 and 95. This means that within their parameterization, the ECC concept has been implicitly taken into account only to some degree.⁷³ Therefore, it is legitimate to consider scaling factors for ions lying between 0.75 and 1 depending on the employed water model so that one avoids overscaling due to partial double-counting for the electronic polarization effects both within the ion and water models.^63,71 This wide range of scaling factors questions the compatibility of non-polarizable water models with the ECC framework. Although some water models have been shown to perform better than others within the ECC framework,⁶³ there is not yet enough evidence to draw general conclusions. Similar problems arise when considering charged groups of amino acids, lipids, nucleic acid, and other readily parameterized components of biological systems. Depending on the particular force field parameterization, electronic polarization effects are already implicitly included to varying degrees and one should again avoid overscaling. In our case, the strategy so far has been to make only minimal modifications to the existing biomolecular force fields, focusing the scaling of the partial charges to the charged (or zwitterionic) groups, with a factor of 0.75 or slightly higher.^4–6,10,72

ECC users should be aware that some quantities, such as the solvation free energy often used in force field parameterization or validation, are not so straightforwardly obtained with scaled charge ECC force fields, and care must be taken in their evaluation. Indeed, as pointed out and discussed early on by Leontyev and Stuchebrukhov,⁵³ the solvation free energy ΔG_solv in the ECC framework should be calculated in two steps—the usual nuclear part usually obtained with alchemical free energy simulations, ΔG_MD, and the electronic solvation energy part, ΔG_el, which can be obtained by solving the Poisson equation or, for simple ions, estimated with a Born solvation model. The total solvation free energy is finally obtained as the sum of these two contributions ΔG_solv = ΔG_MD + ΔG_el, which are of comparable magnitude.

The ECC framework has been first applied to simulations of simple salt solutions made of monoatomic ions. Following the original prescription of the ECC theory, the total ionic charge is fixed at 0.75 times its nominal value (i.e., ±0.75 for monovalent ions and ±1.5 for divalent ions) so that only the two parameters of the 12-6 Lennard-Jones potential remain to be determined for each ion, i.e., the depth of the potential energy well, ϵ, and the distance at which the potential is zero, σ. These parameters need to be slightly adjusted compared to full charge force fields because the charge reduction affects, among others, the ion–water interaction and thus the hydration structure. As a general rule, a 5%–10% decrease in the ion radius allows us to recover the correct hydration structure; however, the exact value is water model dependent. Because of the charge neutrality of salt solutions, it is impossible to obtain experimental data for a single cation or anion. To overcome this difficulty, one strategy is to start with simple salts that are known to exhibit very little pairing in solution (at least at not too high concentrations) and to use local experimental probes to parameterize the ions. In this spirit, force fields for simple monovalent cations and anions such as Na⁺,⁷⁴ Li⁺,⁷⁵ and Cl⁻⁷⁶ were developed using data on the local ion hydration structure obtained from neutron scattering with isotopic substitution.⁷⁷ This strategy improved the ion pairing behavior in such ionic solutions and was also shown to successfully capture the affinity of halide anions to a water–oil interface.⁶⁴ 

While, in principle, the charge scaling factor is fixed in the ECC theory, several groups considered it as an additional fitting parameter,^71,78,79 which may reflect the limitations of water force fields.⁶³ This resulted in suggested scaling factors lying in the 0.85–0.94 interval, depending on the employed water force field and experimental observables targeted within the fitting procedure.^63,71,79,80 As a result, several scaled charge force fields for NaCl in solution are now available.^74,78,80

Scaled-charge force fields were shown to provide a faithful description of simple electrolytes over a broad range of concentrations. Moreover, compared to standard non-polarizable force fields, they significantly improved the description of a number of structural and dynamic properties. For instance, the Madrid-2019 scaled charge force field for a series of simple monovalent ions was shown to properly capture a number of thermodynamic (activity coefficients and solubility) as well as dynamic (viscosity and diffusion coefficients) properties.⁷¹ It also improved the description of the salting out effect in NaCl solutions.⁸¹ In general, scaled charge force fields have significantly improved the description of ion pairing, as well as diffusion, viscosity of electrolyte solutions,^63,71 and water reorientation dynamics in the ionic hydration shell.⁸² 

In the light of these promising results, the ECC scaled charge strategy was then extended toward the more challenging cases of non-transition metal divalent cations that exhibit very strong binding artifacts in simulations. The first cation under study was Ca²⁺, which is highly relevant for many biological processes. The ECC calcium force field was parameterized against experimental neutron scattering data on concentrated calcium chloride solutions that are very sensitive to calcium–water distances.⁸³ The ECC force field was shown to properly describe the structure, viscosity, and diffusion of calcium chloride solutions. One should note here that the use of scaled charges does not always reduce the amount of ion pairing, which is determined by the balance between ion–ion and ion–water interactions. In the case of CaCl₂, the use of ECC force fields indeed leads to a slightly reduced amount of ion pairs in solution. Our ECC calcium force field was later further refined to yield better agreement with both experiments—Extended X-Ray Absorption Fine Structure (EXAFS) and neutron scattering—and benchmark computations (ab initio molecular dynamics).⁸⁴ This approach has also led to a scaled charge Mg²⁺ force field,⁸⁵ which was thoroughly checked against other Mg²⁺ models recently.⁸⁶ A similar strategy was also attempted for Zn²⁺.⁸⁵ Even though Zn²⁺ is not formally considered to be a transition metal, it does form partially covalent interactions with solvent molecules, counterions, and other ligands, which are beyond what ECC (or even explicitly polarizable polarizable force fields) is able to describe. Attempts were made to treat such cases using specific interaction-site potentials.⁸⁷ 

Scaled charge force fields were also recently proposed⁸⁸ for a wide range of polyvalent metal cations, including high valence transition metals such as Y³⁺, Zr⁴⁺, and Hf⁴⁺, using the hydration free-energy and ion–water oxygen distance as target properties. Even if such quantities can be reasonably well reproduced, one should be cautious when modeling highly charged ions with classical force fields, as their specific chemical interactions may not be adequately captured.

For full charge models of ions, Li and Merz demonstrated that Lennard-Jones parameters optimized for the ions within a standard 12-6 interaction potential were able to reproduce either the solvation free energies or the ion–water oxygen distances.³⁹ With the extension to 12-6-4 potential form, these two measurables could be reproduced simultaneously, with an open question concerning the description of ion pairing. In this context, it is worth asking whether ECC, which was developed to reproduce the structural properties of electrolyte solutions, also describes ion–water energetics accurately. To this end, we evaluated (Fig. 1) for some of our cation ECC force fields the solvation free energy ΔG_solv, following the ECC prescriptions (see Sec. II C). The two contributions to the solvation free energy, ΔG_MD and ΔG_el, are comparable in magnitude. The total ΔG_solv was found to be in excellent agreement with experimental data, with similar accuracy as standard force fields, even though it never entered in the parameterization procedure. While such calculations represent a useful validation of the force field, one should keep in mind that the value of ΔG_el, which accounts for half the total free energy, is highly sensitive to the details of the cavity used for its estimation and, thus, has a large uncertainty. Hence, the total solvation free energy could hardly be employed as a target property for parameterization.

The description of salt solutions at the air–water interface is problematic within the ECC framework due to the discontinuity in the high-frequency dielectric constant, as discussed above. As a result, the experimental dependence of surface tension on the electrolyte concentration is not properly captured using ECC ions (see Fig. 2). Here, the surface tension of both tested water models, SPC/E and four-point optimal point charge (OPC), shows too weak dependence on the solution molality. The vertical shift between the models stems from the ability of OPC to better describe the surface tension of the pure water–air interface. This limitation of ECC in describing interfacial behavior has to be taken into account when modeling phenomena involving the air–water interface. It is fortunate that such interfaces are rare in biological systems.

Since ECC prescribes the total charge of the ions in simulations with implicit electronic parameterization, only the two parameters of the Lennard-Jones potential remain to be optimized for monoatomic ions. The situation is significantly more complex for polyatomic ions, as there are many more parameters to determine. To a first approximation, all bonded terms (bonds, angles, and dihedrals) can be adapted without change from regular force fields. In some cases, however, the dihedral terms may have to be adjusted as scaling the charges affects the interactions of atoms separated by three bonds (i.e., the 1–4 interactions). Only the total ion charge is prescribed by ECC, which leaves the partial charges and Lennard-Jones parameters of the individual atoms to be determined.

Many common force fields, including those of the Amber family, use partial charges that are determined by fitting the electrostatic potential obtained by quantum chemical calculations under certain restraints. This provides a charge distribution that reproduces the electrostatic potential around the molecule of interest. The idea behind ECC (see above and the recent review⁶²) is that the electronic polarization effectively reduces this electrostatic potential. Hence, in a most straightforward way, ECC partial charges may be obtained by simply scaling the charges fitted to the electrostatic potential by a factor of $1/ϵel≃0.75$⁠.

The Lennard-Jones parameters still remain to be determined. For monoatomic ions, this has been done such as to reproduce the location of the first peak of the ion–water oxygen radial distribution function, obtained from neutron scattering or EXAFS. For complex ions, such experimental radial distribution functions are rarely known. While standard force fields often use the solvation free energy as a target property, this is not practical for scaled charge models, as discussed above. Alternative targets can be used, such as osmotic coefficients or radial distribution functions calculated from ab initio simulations, but there is a priori no unique recipe. In a first approximation and in view of the scarcity of the available data, many scaled charge force fields for complex ions simply did not modify the Lennard-Jones parameters.

Adopting this strategy for carbonate, nitrate,⁹³ and sulfate⁹⁹ allowed for capturing successfully the structure of concentrated KNO₃, K₂CO₃ and Na₂SO₄ solutions in which ions exhibit excessive pairing and aggregation when described with non-polarizable force fields. A similar strategy was more recently adopted also for the acetate anion.⁸⁴ In this case, it turned out that already after charge scaling, the radial distribution functions between the acetate oxygen atoms and water were in good agreement with those obtained from ab initio molecular dynamics simulations. Thus, further scaling of the Lennard-Jones parameter σ was not required. In addition, this scaled charge acetate combined with the scaled charge calcium developed previously was shown to describe well the structure of concentrated calcium acetate solutions as validated against neutron scattering measurements.⁸⁴ Force fields with either full charges (or with only the cation charges scaled) strongly overestimate ion pairing and aggregation in the solution.⁸⁴ However, it is worth noting that a careful comparison with various experimental data suggests that, while full charge force fields indeed strongly overestimate the binding, the ECC force field with a scaling factor of 0.75 appears to underestimate the binding strength.¹⁰⁰ Force fields based on milder scaling factors, i.e., 0.80–0.85, seem to provide a better agreement with experimental data.^63,71

The same ECC acetate was used together with a scaled imidazolium force field to investigate ion pairing in imidazolium–acetate solutions.¹⁰¹ These two moieties represent a simple model of interacting groups involved in salt bridges in proteins—acetate is a proxy for carboxylate groups in aspartate or glutamate residues and the imidazolium cation is a proxy for the arginine side chain (Fig. 3). While the original full charge force fields overestimated ion pairing, charge scaling led to a structure of the solution in good agreement with the results of neutron diffraction experiments with double isotopic substitution, which provide a sensitive measure of the degree of ion pairing in the solution. The reliability of the ECC force field thus being validated, additional information about the geometry of the ion pair was obtained, such as the existence of an out-of-plane binding mode involving the hydrophobic group of acetate in addition to the standard in-plane geometry. Analysis of the Protein Data Bank revealed that such an out-of-plane interaction motif, although less common than the in-plane one, is found in 25%–55% of protein structures. This suggests that despite the apparent simplicity of such model salt solutions, valuable insights can be gained into the strength and geometrical aspects of biologically relevant interactions, thanks to the use of appropriately scaled force fields. Furthermore, this study confirmed that the strength of salt bridges is strongly overestimated in standard non-polarizable force fields.^15,16

Other polyatomic biologically relevant ions such as multivalent phosphates¹⁰² or hydrogenoxalate and oxalate¹⁰³ anions were also recently successfully improved using the ECC framework.

Many biomolecules such as lipids, proteins, or DNA/RNA are either charged or contain charged groups. The ECC approach for monoatomic and polyatomic aqueous solutions can be extended to more complex solutions and biomaterials¹⁷ as the variation in ϵ_el is small (see Table I). As a result, the charge scaling factor, although strictly speaking environment-dependent, can be well approximated as being constant in these biomolecular systems. This approximation allows the usage of a single scaling factor in the development of force fields where molecular interactions generated for diverse molecules are expected to be mutually compatible.

All cells, organelles, and vesicular bodies are encapsulated by membranes. Their basic constituents are amphipathic lipid molecules. Independently of their net charge, most lipid moieties and their derivatives contain positively and negatively charged groups as well as some highly polar groups. These charged groups include, for example, (methylated) ammonium, choline, carboxylate, and phosphate groups. As with simple molecules containing these motifs, lipids display excessive binding when interacting with other charged groups when described with standard non-polarizable force fields. In fact, the interaction of lipids with charged moieties is a conundrum of great relevance.¹⁰⁸ 

Until recently, the consensus was that cations, even the monovalent ones such as Na⁺, interact strongly with the membranes. This idea originated from molecular dynamics simulations with standard non-polarizable force fields, within which most biologically relevant cations accumulated at and often overcharged even membranes made of only zwitterionic lipids with zero net charge.¹⁰⁹ Yet, the conformation of the lipid headgroup region, as characterized by nuclear magnetic resonance (NMR)-derived order parameters, showed no major changes during titration experiments with monovalent cations.¹¹⁰ A more recent study employing electrophoresis and isothermal titration calorimetry concluded that these simulations exaggerated Na⁺ adsorption onto bilayers composed of zwitterionic phosphatidylcholine lipids.¹¹¹ For divalent ions, the situation is even more dramatic—in standard simulations, Ca²⁺ ions bind irreversibly to membranes consisting of phosphatidylcholine lipids and their strong accumulation at the membrane surface leads to overcharging, which contradicts with most of the available experimental data.⁹ This excessive binding of cations to membranes in all standard non-polarizable force field simulations has been thoroughly addressed recently.^32,112

To address the overbinding problem, the developers of the CHARMM36 force field introduced a specific pairwise non-bonded repulsive potential between certain charged groups to prevent their close approach.³² While this potential results in somewhat extended distances between these groups, it leads to a better agreement with several other experimental measurables. It effectively reduces the interaction and helps to address the overbinding issue. However, this approach, which lacks a clear physical justification, has its drawbacks. Processes that involve polycoordination of ions with several lipids are not described properly. This is especially important for divalent cations such as Ca²⁺ that can mediate binding between like charge regions of signaling proteins and membrane lipids. An example of such a case is the membrane-binding C2 domain of protein kinase C α (PKCα-C2), which is bridged by Ca²⁺ to a negatively charged phosphatidylserine lipid (see Fig. 4). Modeling such situations is beyond the capabilities of the original CHARMM36 force field and its NBFIX modification. As shown in the upper left panel of Fig. 4, a membrane consisting of phosphatidylcholine and phosphatidylserine lipids gets overcharged by Ca²⁺, and PKCα-C2 does not attach to the membrane when the CHARMM36 force field is used without NBFIX. When NBFIX terms are included, PKCα-C2 binds to the membrane. However, the binding is not bridged by Ca²⁺, and the protein orientation does not agree with the one in the crystal structure.

The source of the excessive interaction between ions and lipids in classical force fields is again the lack of electronic polarizability and the associated shielding. Therefore, the use of polarizable force fields should, in principle, solve the problem. While there is a significant progress in the accuracy of polarizable lipid force fields,^21,114–117 their results are not always in better agreement with experiment than those of classical force fields.²¹ Combined with the limited availability of parameters for lipid moieties and a tedious parameterization process, the larger requirements of computational resources render polarizable simulations of limited use in studies of complex biological processes. It is, therefore, tempting to instead correct the classical force fields within the ECC approach.

The very first, to the best of our knowledge, atomistic simulation of a lipid bilayer⁵⁰ used charges reduced by 50% based on a pioneering work on surfactants by Jönsson, Edholm, and Teleman.⁴⁹ However, the scaling factor of 0.5 lacks a firm physical basis. The first attempts to solve the overbinding problem within the ECC framework used scaled ions^74,83 together with unmodified traditional lipid force fields.^8,9 From these studies, three main conclusions can be extracted. First, scaling only the ions is not sufficient to avoid the overbinding of cations to a phosphatidylcholine (PC) membrane.⁹ Second, the change in binding brought about by ECC ions depends on the lipid force field. As an example, ion scaling leads to a decrease in Na⁺ binding to POPC membranes with both Slipids¹¹⁸ and CHARMM36 force fields.¹¹⁹ However, for Ca²⁺, the expected decrease in binding with ECC ions is only observed with CHARMM36, where the binding is no longer irreversible at the simulation timescale. In contrast, the binding is still irreversible and (counterintuitively) enhanced when the Slipids force field¹²⁰ is used. Third, anions such as Cl⁻ do not bind to POPC membranes.

Solving the overbinding problem within the ECC framework requires simultaneous scaling of ions and lipids. This includes the ubiquitous phosphatidylcholine lipids that contain zwitterionic headgroups but have a zero net charge. The ECC framework does not provide a recipe on how the scaling of such molecules should be done. One can scale only the charged and polar groups within the headgroup region, i.e., amine, choline, carboxyl, phosphate, and carbonyl groups, leaving any other apolar groups of the lipid untouched. Alternatively, the charges of the whole headgroup and the glycerol region could be scaled. Both approaches will ultimately decrease the charge separation and local dipoles within a lipid. The best scaling strategy depends on the distribution of partial charges in the original lipid force field.

For lipid models of the Amber force field family, the net charges of charged groups are distributed all over the headgroup and glycerol region. This leads to a significantly smaller charge separation in the headgroup compared to CHARMM-based force fields. To some extent, such charge delocalization resembles the ECC scaling process and may explain why applying the ECC framework is not straightforward for the Amber force field family. Notably, one needs to avoid overscaling of the already (implicitly) scaled charges.

A uniform scaling of the entire headgroup and glycerol regions with a factor of 0.80 was used for the Lipid14 palmitoyloleoylphosphatidylcholine (POPC) model, which led to a dramatic improvement of the binding behavior of monovalent Na⁺ and divalent Ca²⁺ cations.¹⁰ This ECC-POPC model achieves a nearly quantitative agreement with NMR experiments. To achieve such a level of accuracy, the scaled atoms had their Lennard-Jones σ parameters scaled by a factor of 0.89. The same scaling strategy was applied to palmitoyloleoylphosphatidylserine (POPS) using the Amber Lipid17 force field as a template.⁷² In contrast to POPC, POPS is a charged lipid, and thus, the total charge of the ECC-POPS model was scaled to −0.75. Again, the ECC scaling rendered the interactions between POPS and the cations more realistic than in the original force field.

The strategy of scaling the headgroup and the glycerol as a whole was also applied to the CHARMM36 POPC force field, resulting in an improvement of the binding of monovalent and divalent ions.¹⁰ Note, however, that lipid models using the variants of the CHARMM force fields are modular, meaning that a molecule is composed of fragments usually involving functional groups with a net integer charge. In this case, a scaling strategy where only charged fragments are scaled seems to be more appropriate. Such an approach avoids undesired structural changes and minimizes the changes in the interactions with other membrane moieties with respect to the original model. Within this strategy, the nominal charge of each fragment after scaling is the original charge of the group times the ECC scaling factor. For example, a entire choline group with a total charge of +1 e has a total charge of +0.75 e after scaling. Yet, the partial charges of independent atoms do not need to follow this relation. In addition, some strong dipoles in very polar functional groups may need to be scaled, such as the carbonyls resulting from the esterification of the glycerol group and the acyl chains. This is usually done to keep the partial charges on these carbonyl oxygens below the values of their more polar counterparts such as phosphates and carboxylic groups after their scaling.¹⁰ 

Overall, adapting the existing lipid force fields to the ECC framework dramatically improves the description of lipid–ion interaction.^9,10,72 This is achieved without compromising the excellent performance of the original force fields concerning the structural properties of the lipid bilayers. By limiting the ECC scaling only on charged (and possibly also highly polar) groups, accurate binding of lipids to charged moieties is achieved without the need to reparameterize the entire lipid model. The remaining question to address in the future is whether charge scaling also improves the binding of other charged biomolecules, such as proteins, to membranes.

Five out of the 20 essential amino acids—Glu, Asp, Lys, Arg, and His (when doubly protonated)—have charged side chains and play a critical role in protein structure and function. Charged amino acids drive protein folding by promoting their exposure to the solvent. For example, protein stability and structure is dependent on pH, which changes the protonation state of these amino acids.¹²¹ Salt bridges, where a positive and a negative amino acid directly interact, constitute fundamental intra- and inter-protein interactions.¹²² They drive multiple phenomena, such as multimerization and the formation of complexes required to form functional assemblies. Charged amino acids are also fundamental in the binding of ligands and cofactors where they foster the recognition, selectivity, and strength of the binding.¹²³ Furthermore, both monoatomic and polyatomic ions act themselves as protein ligands and cofactors.¹²⁴ 

Considering that the high-frequency dielectric constant (ϵ_el) of the protein environment is comparable to that of water (see Table I), it is not surprising that charge–charge interactions involving proteins are also overestimated when using standard non-polarizable molecular dynamics simulations. This issue was noted early on by people simulating ion channels, where the ions, especially divalent cations, usually get trapped in the protein’s selectivity filter.¹²⁵ For this reason, charge scaling of the ions and even the charged groups in the protein’s selectivity filter is one of the strategies proposed to overcome such excessive interaction and to allow for realistic movement of ions through the channel.⁵¹ The potential of applying charge scaling to amino acids to improve overbinding problems was soon realized by Stuchebrukhov et al.^61,126 when developing the ECC framework. In particular, the description of the important salt bridges in proteins improved dramatically when using the ECC correction.^61,101 ECC-scaled proteins were then adopted and proved extremely useful in studies on systems where the proper description of ion–protein interaction is central.^4–6,127

In particular, Ca²⁺ plays a key role in cellular signaling, one of its most important receptors being calmodulin. Calmodulin has four Ca²⁺ binding sites in its EF hands, each with a different but reversible binding affinity. Standard force fields lead to grossly overestimated binding energies of up to ∼−19 kcal/mol, whereas experiments predict ∼−5.4 kcal/mol.⁴ Using an ECC Ca²⁺ ion model together with the original Amber protein force field decreased the largest binding free energy to −10 kcal/mol, which is still twice the experimental value. Only when the negatively charged side chains were scaled, too, did the binding free energies drop close to the experimental values. Quantitatively, scaling has varying effects on each of the binding site. The binding behavior of Ca²⁺ to calmodulin in a aqueous solution with 50 mM CaCl₂ with and without ECC scaling is illustrated in Fig. 5. It is of importance to highlight the overbinding of Ca²⁺ to calmodulin when the standard CHARMM36 force field without NBFIX is used. We observed that ∼20 Ca²⁺ were bound to calmodulin instead of the expected 4 (which corresponds to one per binding site). We saw that certain calmodulin binding loops embraced two Ca²⁺ ions. Furthermore, even a single carboxylate group occasionally bound two Ca²⁺. This is clearly an artifact, as carboxylate groups are abundant in water-soluble proteins and not all of them bind Ca²⁺. When using the same CHARMM36 force field with ECC parameters, the Ca²⁺ binding is reduced drastically to an average of 6 Ca²⁺ bound to calmodulin. This improvement was not without collateral problems, as the binding affinity to some loops became too low so that they unloaded during the simulations (see the black loop on the left snapshot of Fig. 5). Despite this drawback, ECC rendered the binding profile more realistic, as the binding to carboxylate groups was reduced and a complete Ca²⁺-binding loop was thus required for the interaction to take place.

For another Ca²⁺ binding protein, recoverin, ECC scaling not only improved the binding of Ca²⁺ but also captured the overall protein conformational change upon ion binding. This conformational change is of biological interest as it is required for the activation of the myristoyl switch.¹²⁷ 

Intrinsically disordered proteins, especially highly charged ones, constitute another group of proteins for which ECC scaling provides a promising modeling approach. Most of the commonly used protein force fields predict too stiff protein structures with little room for conformational changes. This results from the traditional need to preserve the crystallographic secondary and tertiary structures of proteins. While many factors may contribute to such stiffness, the interaction between neighboring charged groups is an important factor in some common motifs. This interaction is, in turn, strongly dependent on the ionic environment. For example, the structures of polyglutamic and polyaspartic peptides in solution change drastically in molecular dynamics simulations with the underlying force field and the used counterions.⁶ While the CHARMM27 and Amber ff99SB-ILDN force fields resulted in collapsed structures, the CHARMM36 force field yielded a more extended configuration. The OPLS force field provided an intermediate situation. The use of NBFIX always promoted more extended peptide structures by reducing the ions bound to the peptide. When applying ECC corrections instead, it was evident that the results—more extended or collapsed peptides—depended non-trivially on the original force field.⁶ Again, as with the lipids, the outcome of applying ECC scaling to the CHARMM36 force field is fairly predictable, resembling that of NBFIX. It reduces the number of counterion–peptide contacts. As bound ions shield the intra-peptide repulsion, ECC promotes more extended peptide structures. The overall repulsion between neighboring negatively charged groups is increased, although ECC weakens this interaction too. With the Amber force field, which implicitly accounts for electronic polarization to some extent with its delocalized partial charges, the results are system-dependent. In some cases, the binding of counterions to the peptide is even enhanced, resulting in more compact conformations. Therefore, ECC scaling may alter the balance between intra-peptide and peptide–ion interactions in largely unpredictable ways.

As mentioned in Sec. V A, in some cases, ions can act as bridges for binding of proteins to lipids. Such a case is a demanding benchmark for the ECC approach, since it requires the correct balance between lipid–ion, lipid–protein, and ion–protein interactions. As shown in the top right panel in Fig. 4, our ECC simulations of PKCα-C2 capture its binding to a membrane consisting of POPC and POPS lipids in the presence of CaCl₂. As shown in the bottom left panel of Fig. 4, the binding mode agrees well with the crystal structure, with Ca²⁺ bridging the PKCα-C2 to a POPS lipid. The bottom right panel of Fig. 4 shows the histogram of minimum distances between the center of mass of the Ca²⁺-binding loops of PKCα-C2 and membrane POPS. The ECC simulations are in much better agreement with experiment (gray dashed line) than full charge simulations using CHARMM36 with or without NBFIX. Still, the agreement is not perfect, as the orientation of a POPS lipid differs somewhat from the crystal structure, and for the majority of the simulation, only one Ca²⁺ is involved in the bridging instead of the experimental value of 2 (PDB:1DSY). Moreover, the aforementioned PS–loops distance is slightly larger than in the crystal structure. These differences may have arisen from limited sampling in the simulations or from the differences between the binding mode in the crystal structure and in the native environment. Finally, more work needs to be done to deal with other charged groups found in proteins such as those originating from phosphorylation.

Description of other highly charged biomolecules should also benefit from ECC scaling. Prime examples are sugars such as glycosaminoglycans and sialic acids in N-glycans, as well as DNA and RNA. Applying the ECC strategy to such biopolymers shares the problems discussed for polyatomic systems earlier. So far, there is no unique procedure to derive the scaled charges, which ultimately depends on how the unscaled partial charges were originally derived and distributed within the molecule. Nevertheless, promises of the ECC approach to describe phosphate-containing species have been demonstrated in a recent work on solutions containing divalent and trivalent phosphates, osmotic properties of which were greatly improved by the use of scaled charges.¹⁰² The extension of this approach to monovalent phosphates should be straightforward, but the challenge lies in its application to RNA and DNA. Indeed, in contrast to proteins—where the charges are relatively far from the backbone that defines the secondary structure—in DNA and RNA, as well as in sugars, the charges are located directly on the backbone or sit next to the rings. Therefore, any change in partial charges is expected to have a large effect on the dihedral terms. These dihedral terms have been carefully optimized in full charge force fields to reproduce the complex conformational phase space of these biomolecules. Hence, extension of ECC to these biomolecules without compromising the description of their structural parameters remains challenging.

Despite the apparent simplicity of the ECC theory, which prescribes that partial charges of the charged interacting moieties have to be adequately scaled, its practical implementation is not always that straightforward. First, an increasing number of studies suggest that the theoretical ECC scaling factor of 0.75 often results in too weak charge–charge interactions. Several groups have proposed scaling factors closer to 1 to resolve this issue,⁷¹ which can, in part, be justified by the use of non-optimal water models within the ECC framework.⁶² Second, scaling of charges affects the interactions of the charged group with other moieties such as water, which were often fitted to reproduce experimental data (for instance, ion–water oxygen radial distribution functions). Hence, even for simple salt solutions, charge scaling often needs to be accompanied with refining the Lennard-Jones parameters to restore full agreement with experiment in these properties.¹⁰ More generally, a concrete implementation of the charge scaling depends on how the starting underlying full charges were originally developed for polyatomic groups. Note that there is not a unique set of charges compatible with the ECC framework in the same way as different full charge force fields adopt different distributions of partial charges. From our experience, the charge distributions in several of the majorly used biomolecular force fields already implicitly include the effect of electronic polarizability to a certain extent. In such cases, one should thus be careful not to overscale the charges. The water model, and especially its dielectric constant, has a non-negligible effect on the ECC parameterization, with the same being also true for full charge models, as illustrated by the fact that many recent ion force fields are available, in particular, flavors adapted to each water model. In practice, simply scaling the charges of an existing biomolecular force field is often not sufficient to match experimental data for ion binding.

In addition, one should note that the outcome of charge scaling within the ECC framework, in contrast to some other schemes developed to quench specific charge–charge interactions, is not always the weakening of ion pairing. In particular, because ion pairing is a property that depends on the balance between ion–ion and ion–water interactions, scaling down the ionic charges can sometimes (rarely in fact) result in enhanced ion pairing. In addition, charge scaling may require introducing correction terms for pressure calculation.¹²⁸ Finally, certain limitations of the ECC approach are inherent to its mean-field nature. For instance, ECC can hardly capture cross polarization between multiple charged groups, which may be important in specific binding environments.²⁹ The mean-field nature of ECC also results in inherent issues when dealing with interfaces separating regions with different ϵ_el (see Fig. 2) or when simulating a system immersed in an external electric field.^57,58

Force field molecular dynamics is an efficient simulation technique that has proven its value in providing understanding of molecular mechanisms. Nevertheless, there is increasing amount of evidence suggesting that standard non-polarizable force fields overestimate interactions between charged groups, which are of paramount importance for many biological systems. The central problem is the lack of effective charge shielding due to the electronic polarizability induced by the surrounding environment. One solution to this problem is to use explicitly polarizable force fields. However, issues connected with parameterizations as well as supra-microsecond timescales and large system sizes inherent in many biological problems represent a challenge for polarizable models due to their computational cost.

Several approaches have been developed to attenuate the binding between charged moieties and groups within the classical force field framework, such as tailoring specific pairwise interactions between charged moieties or including electronic polarization in a mean-field way within the ECC framework. The latter, which has been reviewed here, is a physically well-justified approach. In its charge scaling variant, ECC is least intrusive, as only charges on atoms or groups contributing significantly to charge–charge interactions are scaled. This means that partial charges of weakly polarized groups can be safely left unmodified. This approach can be adopted with only minor modifications to the existing biomolecular force fields and thus benefits from decades of their developments. We have shown here that, despite all current limitations, this approach has been very successful, allowing not only for a faithful description of electrolyte solutions but also improving significantly the behavior of charged species in biosystems. Specifically, it leads to the correct description of ion binding to proteins and salt bridges as well as to significant improvements in lipid–ion and lipid–protein interactions, which demonstrates the potential of the ECC approach for a realistic description of complex biomolecular systems.

At present, there is no unique set of charges compatible with the ECC framework. Therefore, additional and non-standardized criteria are required. This calls for a unified protocol for force field development within the ECC framework, possibly going back to the drawing board to develop a de novo consistent and unified model of water, ions, peptides and proteins, lipids, sugars, and nucleic acids. If we find enough stamina for this immense task, this should be the focus of future development of the ECC “ecosystem.”

The simulation data that support the new findings of this Perspective are openly available in the Zenodo repository as datasets at the following DOIs: 10.5281/zenodo.3888488 (Solvation free energy of cations), 10.5281/zenodo.3888369 (Surface tension of ECC NaCl and CaCl₂ models with OPC water model), 10.5281/zenodo.3888383 (Surface tension of ECC NaCl and CaCl₂ models with SPC/E water model), 10.5281/zenodo.3888436 (Surface tension of Dang NaCl and CaCl₂ models with SPC/E water model), 10.5281/zenodo.3891388 (PKCα-C2 interacting with membrane), and 10.5281/zenodo.3888332 (Calmodulin in water with a large CaCl₂ concentration).

E.D.-D. acknowledges support from the “Initiative d’Excellence” program from the French State (Grant Nos. “DYNAMO,” ANR-11-LABX-0011, and “CACSICE,” ANR-11-EQPX-0008). H.M.-S. acknowledges the Czech Science Foundation (Grant No. 19-19561S). P.J. acknowledges support from the Czech Science Foundation (EXPRO Grant No. 19-26854X). M.J. thanks the Emil Aaltonen Foundation for funding. P.D. acknowledges support and computational resources from the European Regional Development Fund OP RDE (Project ChemBioDrug No. CZ.02.1.01/0.0/0.0/16_019/0000729). We thank the IT4Innovations National Supercomputing Center in Ostrava (Project No. OPEN-8-1) and CSC—IT Center for Science (Espoo, Finland) for computational resources.