PCCP Dynamic Article Links A nline Cite this: Phys. Chem. Chem. Phys., 201 1,13,2613-2626 www.rsc.org/pccp PERSPECTIVE Accounting for electronic polarization in non-polarizable force fields Igor Leontyev and Alexei Stuchebrukhov* Received 29th September 2010, Accepted 19th November 2010 DOI: 10.1039/c0cp01971b a. U The issues of electronic polarizability in molecular dynamics simulations are discussed. We argue that the charges of ionized groups in proteins, and charges of ions in conventional non-polarizable force fields such as CHARMM, AMBER, GROMOS, etc should be scaled by a factor about 0.7. Our model explains why a neglect of electronic solvation energy, which typically amounts to about a half of total solvation energy, in non-polarizable simulations with un-scaled charges can produce a correct result; however, the correct solvation energy of ions does not guarantee the correctness of ion-ion pair interactions in many non-polarizable simulations. The inclusion of electronic screening for charged moieties is shown to result in significant changes in protein dynamics and can give rise to new qualitative results compared with the traditional non-polarizable force field simulations. The model also explains the striking difference between the value of water dipole \i ~ 3D reported in recent ab initio and experimental studies with the value \fK ~ 2.3D typically used in the empirical potentials, such as TIP3P or SPC/E. It is shown that the effective dipole of water can be understood as a scaled value \ieS = n/^/Sei, where £el = 1.78 is the electronic (high-frequency) dielectric constant of water. This simple theoretical framework provides important insights into the Department of Chemistry, University of California Davis, One Shields Avenue, Davis, California 95616, USA Igor V. Leontyev received his PhD in 2004 under the mentorship of Prof Mikhail Basilevsky from Karpov Institute of Physical Chemistry (Moscow, Russia). PhD thesis: "Molecular mechanisms of equilibrium and non-equilibrium solvation effects." After completion of his PhD degree he spent a year in the laboratory of Prof Masanory Tachiya at National Institute of Advanced Industrial Science and Technology (Tsukuba, Japan) to study intermolecular electron transfer reactions by polarizable molecular dynamics. In 2006 he joined Alexei Stuchebrukhov's group as a postdoctoral fellow, and is currently studying the mechanism of respiratory enzyme cytochrome c oxidase using MD methods. Since the success of molecular simulations depends on the quality of the force fields used, he is also conducting research into basic principles of biological simulations, focusing on the electronic polarization issues. Igor Leontyev / Alexei Stuchebrukhov received a PhD degree in theoretical physics from Moscow Physical and Technical Institute in 1985, in the theory group of Academician Vitaly L. Ginzburg, and was a research fellow in the Institute of Laser Spectroscopy and the Research Center for Laser Technology of Russian Academy of Sciences until 1990. During this time he worked in a Laboratory directed by V. S. Alexei Stuchebrukhov Letokhov. In 1990 he joined the group of Rudy Marcus at California Institute of Technology, with whom he was a research fellow until 1994. Until this time his research interests were mainly in the area of vibrational dynamics and spectroscopy of polyatomic molecules. Since 1994 he has been at UC Davis, where he is currently a Professor of Chemistry and Biophysics. His current interests are in the general area of chemical kinetics, and biological electron and proton transfer reactions; his main current efforts are focused on understanding molecular aspects of electron transport chain of aerobic cells and biological energy transduction. He has been a recipient of Alfred P. Sloan Foundation Fellow (1996) and Arnold and Mabel Beckman Foundation Young Investigator Award (1997). This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2613 View Online nature of the effective parameters, which is crucial when the computational models of liquid water are used for simulations in different environments, such as proteins, or for interaction with solutes. a. U 1. Introduction At present, the majority of molecular dynamics simulations are performed by using non-polarizable force fields such as AMBER,1 CHARMM,2 GROMOS3 and OPLS.4 In these models, the all-important effects of electronic polarization and screeningt of electrostatic interactions are presumably incorporated in the effective charges and other empirical parameters of the force fields. Despite the drastic simplifications, non-polarizable models have been remarkably successful in modeling many complex molecular systems.6 For example, the properties of liquid water are described quite accurately without introducing electronic polarizability explicitly; likewise, the hydration free energies can be computed quite accurately using non-polarizable simulations.7'8 However, the simulation of polarization effects in low-polar solvents, e.g. ethers,9 and especially in non-polar solvents, e.g. alkanes,10'11 meet serious problems. The non-polarizable models can also significantly underestimate the magnitude of the dielectric response in low-dielectric protein environment12'13 and lipid membranes.14 For example, the dielectric constant of the inner part of cytochrome c was found to be only about 1.5,15 which is lower than pure electronic dielectric constant £ei s 2.0.16 Many other shortcomings of non-polarizable MD simulations have been recently discussed in the literature, see ref. 17 and references therein. The polarizable models aim at resolving the problems mentioned above. Most of such models involve various kinds of coupled polarizable sites,9'10'18'19-22 and the computationally-expensive procedure of achieving self-consistency of such sites at each molecular dynamics time step .J The implementation of such models is yet to be completed; at present, even the simplest classical Drude oscillator model9'10'20-22 is still not readily available for application to many biological systems. As fully polarizable force fields are being developed, there is also a clear need for better understanding the existing non-polarizable models, in particular how accurately they capture the effects of electronic polarization and screening,23'24 and possibly improving them. Given a specially designed (but empirical in nature) procedure of how the partial charges are selected1'2 the charges of neutral residues do reflect, at least approximately, the effects of electronic screening—in a way how for example TIP3P or similar fixed-charge models of water does so. One issue of concern, however, is that the electrostatic interactions of ions are described in standard non-polarizable force fields, such as CHARMM or AMBER, by their original integer charges (e.g. ±1, for Na+ and C7~), i.e. as if these ions were in vacuum, completely disregarding the effect of electronic dielectric (e = £el) screening inherent to the condensed phase medium. The interaction of such bare t Throughout the paper the term "electronic screening" means a reduction of the electric field and electrostatic interactions due to an electronic relaxation of the environment. For the origin of the effect see e.g. ref. 5. \ With the Extended-Lagrangian technique91018,20-22 the computation cost of polarizable simulations can be significantly reduced. charges obviously is overestimated by a factor of about two (the screening factor £ei is about 2 for most of organic media).10 Thus, for example, in simulation of ion channels with conventional non-polarizable force fields, the direct Coulomb interaction of ions (e.g. several K+ ions in the same channel, just a few angstroms apart)25 is probably twice as strong as it should be; similarly, the interactions of ions with water molecules, or with partial atomic charges of a protein are likely to be overestimated as well. The same is true for interaction of charged residues in proteins, such as Arg+ or GliT, partial charges of which carry their original net values ±1. The use of bare charges in non-polarizable simulations would be appropriate for vacuum, but not for condensed phase, where all charges are essentially immersed in the electronic continuum, which weakens their interactions by a factor of about two. A similar problem arises in QM/MM calculations, where one needs to evaluate the electric field of the protein medium to which the QM system is exposed. The use of CHARMM or AMBER charges in such calculations has become standard, and has been adopted in many studies.26 Obviously, the electric potential of charged residues in such calculations should reflect the electronic screening of the medium. In this paper we discuss a principle of uniform charge-scaling based on which one could systematically build a non-polarizable force field for simulations of condensed media. The principle is based on a simple idea of uniform electronic continuum, with an effective dielectric constant e ~ 2, and point charges moving in it. The resulting model, which combines a non-polarizable (fixed-charge) force field for nuclear dynamics (MD) with a phenomenological electronic continuum (EC) is referred to as MDEC (Molecular Dynamics in Electronic Continuum). In some sense, the model is an opposite limit of fully polarizable models involving polarizable point dipoles. Of course, the reality neither involves point dipoles nor the completely uniform electronic continuum, but lies somewhere in between these two limits. In MDEC model the effects of electronic screening are reduced to simple scaling of partial charges. The model is similar but not equivalent to standard non-polarizable force fields.1'2"* An ultimate rigorous implementation of the new concept, of course, would require a consistent re-parameterization of all force field parameters such as bond-length, angle, torsion and van-der-Waals parameters along with the effective partial charges. In this review, however, we examine only a simple scaling of partial charges, which makes non-polarizable force fields such as AMBER and CHARMM to be uniformly consistent with the idea of electronic screening that naturally should improve the quality of these force fields. Of our particular interest is the calculation of solvation effects using MD. In the MDEC the electronic polarization part of the solvation is calculated explicitly from the electronic continuum model, while the nuclear part is obtained with a fixed-charge MD. The two parts need to be combined to obtain the total solvation energy.27'28 We will demonstrate that MDEC model and the Drude oscillator model produce comparable results for dielectric constants of alcohols and 2614 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011 View Online a. U alkanes. It will be then argued that using this model one can rather accurately describe the solvation effects and dielectric properties of non-polar liquids and proteins. Several examples of MDEC calculations and the effects of electronic polarization will be discussed, including the interaction between Na+ ions, which is of interest for ion-channel simulations, and the dynamics of an important salt-bride in Cytochrome c Oxidase. Another important issue related to non-polarizable models concerns water. We show that such models as TIP3P and SPC/E in effect are MDEC models. Our theory explains the striking difference between the value of bulk water dipole Hi ~ 3D reported in recent ab initio and experimental studies and the value jfa ~ 2.3D of TIP3P or SPC/E. We show that the effective dipole of water can be understood as a scaled value neS = where £ei = 1.78 is the electronic (high-frequency) dielectric constant of water. This simple theoretical framework provides important insights into the nature of the effective parameters, which is crucial when the computational models of liquid water are used for simulations in different environments, such as proteins, or for interaction with solutes. 2. Theory. MDEC model The MDEC model was discussed previously in ref. 27-29. Here we summarize the main features of the model essential for the subsequent discussion. The MDEC model considers point charges moving in homogeneous electronic continuum of known dielectric constant £ei. The interactions between charges in such a system are scaled by a factor 1 /fiel. It is instructive to see how this model arises from a microscopic polarizable model as an approximation. A formal analysis of this is given next. 2.1 Screening effect and effective charges Consider a system of polarizable point charges. The energy of such a system is written as follows: W(ru...,rN) 1 N i mi Z j*i r,J J2 diK(i,j)dj - J2 E\r¥i (2.1) where q's are the gas-phase partial atomic charges, and ds are the induced point dipoles located at the positions r's of the corresponding charges. The dipoles are induced by the electric field from other charges and other dipoles. The dipole-dipole interaction is quadratic and is described by the matrix K; the diagonal elements of this matrix are inverse polarizabilities, 1/a, which are assumed to be the same for all charges. In addition to the dipole-dipole interactions, the dipoles also interact with the electric field of other point charges, E'(r,). The field is taken at the position of a given polarizable site (and corresponding charge) rt, and the prime indicates that the electric field does not include the field of the point charge itself. For simplicity, hereafter, the vector notations are omitted while the usual vector nature of the appropriate variables is assumed; thus, e.g. E'(r,)di denotes a scalar (dot) vector product, and djK(iJ)dj stands for the vector tensor vector product. The polarizable dipoles represent electronic polarizability of the ions, and therefore respond to an external field "instantaneously". The external field here is the field of point atomic charges, which is changing together with the position of the nuclei on a much slower time-scale than the electronic response. Thus the polarization dipoles are always at "equilibrium" (i.e. minimizing the total energy) for a given configuration of the nuclei, and the dynamics of the nuclei coordinates r can be described with a Born-Oppenheimer type of effective potential energy W(rx,.. .,rN)\ the dynamic coordinates of the dipoles are not present explicitly in this picture. The equilibrium values of the dipoles can be found by minimizing the energy with respect to the dipole values. Each of the dipoles will have the following equilibrium value: E'{ri)-Y/K{iJ)dj (2.2) where the first term in parenthesis is the electric field of the charges other than qb and the second is the electric field of other dipoles dj at the position of the dipole dt. All equilibrium values of the dipoles depend self-consistently on each other, and on the position of the nuclei, which determine the "external" field to which the dipoles are subjected to. The substitution of the above equilibrium values for dipoles into energy expression gives W(ru...,rN) 1 /V 1 1 ^ qiqj 1 (2.3) The second term of the above equation is the energy of the dipoles, which is the same as electronic polarization energy. The point dipole polarization is now written in terms of the polarized continuum as follows: P\r)dr (2.4) where P'(r) is polarization density, and the integration is over the volume Vt of the ith atom. Since the boundaries between atoms are not well defined, here already the approximate character of the treatment becomes evident. The prime of the polarization density indicates that this is part of total polarization at point r caused by the electric field other than the field of the atom itself. This polarization is proportional to the local external field, as in the usual macroscopic continuum electrostatics (here the electric displacement D is the same as E'): P'( 1 £el - 1 471 E'( fiel (2.5) Assuming now that the external electric field E'(r) does not change significantly within the atomic dimensions (this is the second major approximation), the polarization energy can be written as follows We\ = --YJE'(rl)dl = - fiel fiel E Ea(r) dr (2.6) This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2615 View Online a. U er er which after some additional transformations and assuming spherical shape of the atoms becomes WA = -\\- 1 n n2 (2.7) where R's are the corresponding radii of the atoms. The substitution of the above relations into eqn (2.3) gives for total energy the following: 1/fiel) 2Ri' (2.8) As expected, the first term, which represents interatomic interactions, is scaled by a factor l/fiel. The second term in this expression is the polarization (Born) energy of individual ions, which does not depend on the atomic configuration and, therefore, is not important for interaction and dynamics of the charges; however it becomes critically important in solvation energy calculations. The above expression could have been written from the start for a system of charges in electronic continuum. The above derivation shows that the system of polarizable dipoles can be approximated by an effective non-polarizable model. It follows from the above derivation that the dielectric constant £ei is related to polarizability of the dipoles a and density of the polarizable sites Nx as given by the Clausius-Massotti expression: fiel - 1 4nct.No, 1 - (4n/3)xNx (2.9) Thus, the energy of a system of point polarizable dipoles can be approximated by that of an equivalent continuum model. The numerical quality of the approximation is difficult to evaluate a priori, however, despite the known steps of the derivation that involve approximations. We will check the quality of this approximation by comparing directly the results of calculations using Drude model and an equivalent continuum model later in the paper. From the above treatment it follows that the system of polarizable point charges can be substituted by a system of non-polarizable point charges of scaled values qf* = qi/^/s^, so that the interaction between the scaled charges correctly reproduces the actual interaction q'fq'fjry = q^ijl&efij as if they were in vacuum. The model described above is strictly valid only for ionic mono-atomic liquids.30 However, if one deals with groups of charges, which represent molecules, e.g. water, instead of individual point charges, the results are formally the same as above; only in this case there is no simple relation between partial charges of molecules in vacuum and effective charges in the condensed phase. A good example of such situation is water molecule discussed in the text. It should be noticed too that in a different type of a model, in which an artificial boundary between the molecule and the rest of the electronic polarizable continuum is introduced, for the case of neutral molecules the screening factor may differ from 1/e and depend on a shape of the molecular cavity§, this difference, however, is not significant, see Appendix. § For example, for point dipoles in spherical cavities, the scaling factor for interaction is [3^/e/(2e + l)]2; for electronic dielectric constant of 2 (1.78 for water) the difference may not be significant, considering the uncertainty of the shape of molecular cavity. See details in Appendix. In general, the magnitude of the relative dielectric constant £ depends on which part of the medium relaxation is considered explicitly (as moving charges q,) and which part is described phenomenologically as a polarizable dielectric.31 Since in non-polarizable microscopic models the atomic motions are described explicitly, the screening factor should include only electronic component of the medium polarization, e = se\. The static (i.e. time-independent) dielectric approximation in this case is quite accurate, because on the time scale of nuclear motion the electronic relaxation occurs almost instantaneously, reducing at once all interatomic electrostatic interactions by a factor of £ei. The phenomenological parameter £ei can be theoretically evaluated by Clausius-Massotti eqn (2.9) or directly measured as a high-frequency dielectric permittivity (£ei = n2, where n is a refracting index of the medium). For organic materials typical values of £ei are in the range 1.7-2.210. Thus, the uniform dielectric approximation with £ei = 2 can be a good approximation for many biological systems. The resulting model, which combines a non-polarizable (fixed-charge) force field for nuclear dynamics (MD) and a phenomenological electronic continuum (EC) for the electronic polarization, is referred to as MDEC.29 Summarizing, MDEC model considers charges qt moving in electronic polarizable continuum of known dielectric constant se\, see Fig. la. In the uniform dielectric all electrostatic interactions are scaled by a factor l/£el. Since interactions are quadratic in charges, the effects of electronic dielectric screening can be taken into account implicitly by using scaled partial charges, q'f = qi/In the section "Applications of MDEC model" we show by the computational tests that the simple scaling procedure results in accurate effective interaction between ions of general shape in real solvents. The un-scaled original charges are difficult to specify a priory in general (they are not the same as partial atomic charges of a molecule in vacuum, see ref. 29), unless one deals with ions or ionized groups in proteins, whose un-scaled net charges are known. But effective charges q'f can be found empirically by fitting experimental data32'33 or appropriately scaled ab initio interaction energies.' 2.2 Solvation free energy In MDEC model, when the solvation free energy of a group is considered, the electronic polarization free energy is treated a) ® ® 0 0 b) 0 ®£=ßel 0 s=eel 0 0 8 0 0 8 0 solvent ® ® • 0 Fig. 1 (a) MDEC model for the electrostatic interactions between solute (large crossed circles) and solvent (small crossed circles) charges moving in the electronic continuum of dielectric constant eei; the same electronic dielectric constant eei is assigned for both the solvent and solute regions; (b) MDEC model for estimation of the pure electronic part of the electrostatic solvation free energy AGe\. 2616 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011 View Online a. U explicitly. The solvation free energy consists of the nuclear part AGnuc, evaluated by MD, and the pure electronic polarization part AGei (which corresponds to the last term in eqn (2.8)) evaluated by using the polarizable continuum model34 (i.e. by solving the Poisson equation with corresponding boundary conditions, with dielectric constant e = 1 inside the solute region and e = se\ outside, as shown in Fig. lb): AG = AGn AGel. (2.10) The origin of the last term can be traced in the derivation of the previous sub-section, where the last term in eqn (2.8) corresponds to AGei for spherical ions. The derivation can be extended to the case of molecules of a general shape, in which case the electronic polarization energy for each molecule a takes a familiar reaction field energy form, AGS (2.11) where /£ei- The total MDEC polarization free energy of the medium then is: AG = —AGmd + AGei (2.12) fiel where AGmd is the electrostatic solvation free energy obtained in non-polarizable MD using un-scaled solute charges^} (i.e. the standard approach), and AGei is the pure electronic part of the free energy. A more detailed description of the free energy simulation technique accounting for the electronic polarization can be found in ref. 27-29. Thus, unlike the common empirical models,1-4 in MDEC model the electronic polarization term appears explicitly in the expression for the solvation free energy (2.12). The electronic component constitutes more than half of the solvation free energy for ions and affects the entire non-polarizable concept, including the parameterization strategy. The reason why the conventional non-polarizable force fields1-4 reproduce the hydration free energies of ions quite accurately, completely ignoring the electronic part of the solvation energy will be explained in the Section "Empirical force fields". K The charges of the solvent molecules (such as water) are assumed to be already scaled, which is the case in standard AMBER,1 CHARMM,2 GROMOS3 or OPLS4 MD simulations; the empirical charges of neutral species, in contrast to charged species, are typically correctly reflect the condensed matter nature of the interaction. 2.3 Dielectric constant of the medium The dielectric constant of the medium is often employed in the continuum electrostatic e.g. for solvation free energy evaluation calculations.36 In microscopic calculations, on the other hand, the solvation free energy is obtained directly from MD simulations. The question arises often as to what is the effective dielectric constant of the medium emd that corresponds to a specific microscopic model of the system. The free energy relationships discussed in the previous section allow one to make a connection between the total (static) dielectric constant, e0, which includes both nuclear and electronic polarization effects, and the dielectric constant of non-polarizable MD simulations, emd, which does not explicitly describe pure electronic polarization of the medium. Suppose we consider a spherical ion or a pair of spherical ions; in this case, according to ref. 37, the solvation energies will be proportional to their corresponding Born factors: AG ~ (1 - 1/fio), AGei ~ (1 — l/fiei) and AGmd ~ (1 — I/emd)- From eqn (2.12) we find: £0 — £MX>'£el (2.13) That is the total dielectric constant of the medium e0 is not equivalent to that reproduced by the (non-polarizable) MD simulation, emd\ instead, the relationship between the two is given by the above formula. Although the above arguments are strictly valid only for spherical ions and for the bulk solvent modeled with periodic boundary conditions,38 the relation (2.13) is in fact more general. The above relation can be also obtained using the well-known expression38 for the static dielectric constant: £0 = £el + An Wk^T (M2) (2.14) Here (M") is the mean square fluctuation of the total dipole of the dielectric sample V; kB and T are Boltzmann constant and temperature, respectively. According to the MDEC scaling procedure, the actual dipole moment \i of particles in the bulk is related to the effective moment \fK of these particles in non- polarizable model as \i ■ Edit r; therefore, (M2) = ee]{M2MD); where (M"MD) is the mean square fluctuation of the dipole moment observed in a non-polarizable MD. Thus, eqn (2.13) is obtained from eqn (2.14) by noticing that emd is defined via fluctuation (M2MD) with £ei = 1 in eqn (2.14). 3. Applications of MDEC model 3.1 Ab initio interactions modeled by charge scaling According to MDEC model the electrostatic interactions between ions in the condensed phase should be reduced (scaled) by the factor se\, electronic dielectric constant, with respect to those in gas-phase. To test how well the simple charge scaling procedure reproduces the screening effect for spherical and non-spherical ions in electronic continuum with appropriate boundary conditions we considered39 interaction of several charged species in ab initio calculations. The ab initio treatment captures the effects of electronic polarization of charged species themselves, while the effects of the polarization of the environment, and corresponding screening, are This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2617 View Online a) Na* - Na* a. U Fig. 2 Interaction energies between (a) Na -Na ions; (b) Arg and GlvT amino-acids; (c) GlvT amino-acid and water. Open circles stand for the energies obtained in gas-phase HF/6-31(d) calculation; filled squares are for the same interactions but calculated in dielectric of £ — 2.0. The dashed lines represent interaction energies obtained by the standard CHARMM force field, and using TIP3P water model in (c). The solid lines represent the interaction energies obtained by the CHARMM force field with scaled charges (£d = 2.0), and TIP3P water model in (c). described here phenomenologically, by a continuum with dielectric £ei = 2.0. The interaction energies were calculated using a quantum-mechanical procedure identical to that of the CHARMM parameterization protocol.2 For each structure with a given separation r between moieties, the interaction energy was calculated as the difference between the total supermolecule energy and the sum of the individual monomer energies. The gas-phase interaction energies were calculated with model compounds in vacuum; while the bulk-phase interactions were obtained with the model compounds immersed in the dielectric of £ = 2. The quantum-mechanical calculation in dielectric utilized the PCM34 technique and self-consistent reaction-field procedure implemented in Gaussian03.40 In Fig. 2a, b and c, the ab initio interactions between ions Na+-Na+, Arg+-GlvT and between GlvT ion and water are compared with those modeled by the original and scaled CHARMM force fields. (For amino acids their corresponding model compounds are used.) In all cases, as expected, there is a significant screening effect of the dielectric environment on the interaction energy; the effect as seen, however, can be pretty accurately reproduced by a simple scaling of charges. Notice that the charges are scaled by a factor 1 /given what have been said about TIP3P water model (see Introduction, and the following sub-section on water), in the example of GlvT-H20 pair, only charges of GlvT were scaled. Notice how accurately the scaled CHARMM force field reproduces the results of ab initio calculations. Similar results are expected for other common force fields (such as AMBER,1 GROMOS,3 etc) where charged groups are also treated as having their original vacuum net charges. Thus, the simple charge scaling procedure in standard non-polarizable force fields can account for the effects of electronic screening not only in the interactions between ions but also between ions and water. Although, as seen in Fig. 2a, some additional adjustments of van der Waals parameters might be useful to improve the interactions at shorter distances. 3.2 Water models Many non-polarizable force fields are essentially MDEC models. For example, TIP3P32 or SPC/E33 and similar models of water involve empirical charges that can be considered as scaled charges. TIP3P is particularly interesting in this regard as it is often used in biological simulations, and it serves as a reference for phenomenological parameters assignment of CHARMM.2 It is known that the dipole moment of a water molecule in vacuum is 1.85D; in liquid state, however, the four hydrogen bonds to which each water molecule is exposed on average strongly polarize the molecule and its dipole moment becomes somewhere in the range || of 2.9D to 3.2D.42'43 The significant increase of the dipole from \ia = 1.85D to a value \i x 3D, or even larger, is also supported by the Kirkwood-Onsager model,44 which estimates the enhanced polarization of a molecule due to the reaction field of the polarized environment. Yet, the dipole moment of TIP3P water model is only 2.35D. The specific value of TIP3P dipole moment can be understood as a scaled dipole, so that the dipole-dipole interactions are screened by the electronic continuum by a factor l/se\. Indeed, if each dipole (or all partial charges) is scaled by a factor 1 / ~JĚě\, one could consider interaction of the effective dipoles, neíí = /x/i/ěď ^ 2.35Z) (for water se\ = 1.78), as if they were in vacuum. This appears to be exactly what the fixed-charge water models do. Thus, the charges of TIP3P water model should be understood as scaled charges that reflect the effect of electronic screening. (The scaling factor 1 /fiei for water dipoles is an approximation, and one can argue that a different scaling factor should be more appropriate, however, numerically all reasonable continuum models give about the same result, see Appendix.) The scaled nature of charges of TIP3P water model becomes critically important when the interaction with a solute is considered. For example, if the charge of say Na + || It is recognized that in ab initio simulations of bulk water the water dipole can not be defined unambiguously and depends on the partitioning scheme used;41 as such, its actual value remains a matter of debate. Here we rely upon calculations and the partitioning scheme of ref. 42. 2618 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011 View Online a. U Fig. 3 Dependence of water dipole on the polarity of the environment, as given by the Kirkwood-Onsager model,44 eqn (3.1). The parameters are /i0 = 1.855D, a — 1.47 A3 and R — 1.55 A. ion is assigned to be +1, then it is obviously inconsistent with the charges of water model; as the latter are scaled by a factor of while the charge of the ion is not. Clearly the strength of interaction is overestimated in this case by a missing factor of l/^/Ěď, i.e. about 0.7 (for proteins £el ~ 2). The problem would not arise if the charge of the ion were appropriately scaled. (The reason why seemingly incorrect charge gives reasonable aqueous solvation free energy is explained in the next section.) A detailed discussion of non-polarizable models such as TIP3P and SPC/E is given in our recent paper ref. 45. Here we briefly consider the transferability issue for water models that can be demonstrated using the Kirkwood-Onsager model.44 In this model, the medium is represented by a continuum dielectric of e0 and the solvent molecule is modeled by a point polarizable dipole, placed in a spherical cavity of radius R; the permanent dipole is /x0 and the polarizability is a. In such a model, for the average dipole moment of bulk water \ii one obtains the following expression: ft ft) , 2(e0 - 1) a ' (2e0 + l)i?3 (3.1) To use this model for estimation of the liquid state water dipole, one needs to know the value of the radius R of the molecular cavity. This is obviously a phenomenological parameter, which needs to be fixed in comparison with experimental data or derived from a suitable theoretical model. One reasonable estimate of the radius of molecular sphere R is to assume that 2R is the average distance between the liquid water molecules, a = 2R. The distance between the molecules is related to their number density, Nx, as Na = a~3 = 1/8R3. Employing now the Clausius-Massotti relation between Nx and £ei as given by eqn (2.9) the radius R is expressed as R' n (fiel + 2) 6(£el-l)C (3.2) With the radius R given by eqn (3.2), the equilibrium dipole moment in the liquid phase becomes ft ftj x 6(^-1)2(80-1)' % (fiel + 2) (2fi0 + 1) (3.3) which gives the value 3.0 D for the bulk water (fiei = 1.78, £0 = 78, He = 1.85D). This value is in good agreement with the recent experimental estimations; however, it is slightly different from the empirical TIP3P or SPC/E MDEC value, 2.35Z>v/1.78 = 3.14D. This slight inconsistency can be corrected by the fine tuning of the model parameters in computational parameterization procedure.45 The above model will now be used to examine the transferability of non-polarizable TIP3P or similar potentials for bio-molecular simulations in conditions which are quite different from that in pure liquid state of water. Although (3.1) is a crude model, it nevertheless qualitatively captures the electronic polarization effect induced by the polarizable environment. Taking jt0 = 1.855D for water dipole in vacuum,46 a. = \A1 A3 for polarizability47 and the water cavity radius estimated in ref. 45 as R = 1.55 A, one obtains an estimate of \i in different environments characterized by dielectric constant £0. The dependence of \i on £0 is shown in Fig. 3. In the high-dielectric region, £0 > 20, as seen, the water polarization is almost constant and similar to that of the water molecule in the bulk, £0 = 80. At smaller £0, however, the model indicates a significant dependence of the water dipole moment on the polarity of the environment. As shown in Fig. 3, the dipole moment of a water molecule in the media with £0 < 20 is significantly lower than the value \it of water in the bulk. Thus, in low-dielectric environments, such as proteins or membranes, water should be modeled using potentials different than those of TIP3P or SPC/E (the work to develop such models is in progress in our group at present). 3.3. Empirical force fields In non- polarizable force fields of AMBER,1 CHARMM,2 GROMOS3 or OPLS4 the atomic partial charges of non-charged groups can be understood approximately as "scaled MDEC charges", because, as discussed above, these empirical parameters were chosen in such a way as to reflect the condensed matter nature (including screening) of the interaction. In contrast, the charges of ionized groups remain un-scaled and therefore do not reflect effects of electronic screening, and as such are treated as if they were in vacuum. The non-polarizable TIP3P potential is particularly interesting, as it is often used in biological simulations, and it serves as a reference for phenomenological parameters assignment of CHARMM.2 The scaled nature of charges of TIP3P water model is important to bear in mind when the interaction of such water models with a solute is considered. For example, if the net charge of say GlvT ionized side chain is assigned to be — 1 in simulation, then it is obviously inconsistent with the charges of water model, as the latter are scaled by a factor of 1 /'v/fieT, while the charges of the ion are not. Clearly the strength of interaction is overestimated in this case by a missing factor 1 /^/%\, i.e. about 0.7 (for proteins £ei ~ 2). The problem would not arise if the charge of the ion were appropriately scaled. In free energy simulations with non-polarizable force fields (and un-scaled solute charges), the pure electronic contribution to the electrostatic free energy is typically completely This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2619 View Online a. U ignored, as e.g. in ref. 7 and 35. Yet, in many cases such simulations pretty accurately reproduce experimental solvation energies; this may appear surprising, given the fact that about a half of the total solvation free energy (for charged solutes typically 50-100 kcal mol-1) comes from electronic polarization of the medium. In fact, the neglect of large electronic polarization free energy is almost completely compensated by the use of "incorrect" bare solute charges in such simulations. This fortuitous compensation of errors, however, occurs only in the high-dielectric media, as can be seen from the following argument. Consider for example Born solvation energy of Na+ ion, Q = +1, in water; in standard simulations one would have approximately 2R \ emd) (3.4) where emd is the dielectric constant of water that corresponds to a specific MD model employed in the calculation. No matter which model of water is used, emd is much larger than unity, hence the overall estimate of the solvation free energy is Q~/2R, which is independent of properties of the solvent, and can match pretty well the experimental value, provided the ionic radius R is chosen correctly. The interaction between two charges consequently is taken to be then Q"/r, completely disregarding the electronic screening of the interaction. MDEC model suggests instead that in MD simulations the charge Q should be scaled, and the electronic solvation free Q2(. 1 energy AGd 2R V 1 £el added explicitly. In this case, the nuclear part of the free energy calculated in MD will be "2/ __u -IR— given by AG = emdJ and the total free energy, eqn (2.10), is 2R 1 -- —) P. 2R 1 -■ fiel (3.5) Since emd = fio/fiei, as given by eqn (2.13), the above expression correctly reproduces the expected result (Q2/2R)(\ — l/fi0). Notice that the charge is not scaled when the solvation is calculated in electronic continuum. Notice also, that it is only when emd » 1, the two eqn (3.4) and (3.5) approximately give the same result. Yet, for interaction energy of two charges the MDEC gives the correct expression Q~/rse\, while the standard approach gives Q~/r. It is seen that in high dielectric medium, the un-scaled charges result in twice as large contribution for the nuclear part of solvation energy compared with the correct value; as a result the missing electronic part, which is half of the total, exactly compensated. Given that the un-scaled relation is only formally correct when EMD ~ fio/fiel » 1> (3.6) it is not surprising that the traditional non-polarizable approach works well in aqueous solutions (fio/fiei ~ 40), as e.g. in ref. 7 and 28; however, the approach fails (i.e. significantly underestimates the polarization effects) in low dielectric media (fio/fiei ~ 1) as in ref. 9,10,13,14,48,49. Table 1 Dielectric constant of bulk alcohols simulated by different MD models at T — 298.15 K Alcohol £0, exp emd, npol MD* £0, pol. MD' Eei, P0l- MD' £0: MDEC MeOH EtOH 2-PrOH 2-BuOH 1-PrOH 1-BuOH 32.61 24.85 19.26 15.94 20.52 17.33 17.2 18.8 13.7 7.8 15.2 10.8 30.1 21.4 17.6 15.8 19.5 21.2 25.8 30.08 23.29 13.26 24.32 18.36 Rmsd of 6e (%) 3.7 0.6 ' Experimental values50; b Conventional non-polarizable MD model;21 c 0.9 Polarizable classical Drude oscillator model; MDEC model, eqn (2.13), where emd and £d are taken from and c, respectively; € The relative error of Born factor, eqn (3.7). Table 2 Dielectric constant of bulk Alkanes simulated by different MD models Alkane T, K £o, exp" £MO, npol MD* £0, pol. MD1 Eel, P0l. MD' £0, MDECrf Ethane 184.55 1.7595 1.014 1.707 1.697 1.721 Propane 231.08 1.7957 1.015 1.798 1.768 1.795 Butane 272.65 1.8098 1.016 1.801 1.774 1.802 Isobutane 261.43 1.8176 1.015 1.905 1.823 1.850 Heptane 298.15 1.9113 1.018 2.021 1.977 2.013 Heptane 312.15 1.8904 1.018 1.976 1.933 1.967 Decane 298.15 1.9846 1.020 2.118 2.066 2.106 Decane 312.15 1.9668 1.019 2.128 2.074 2.113 Rmsd of 6e (%) — — 96.4 5.1 — 4.3 ° Experimental values50; b Conventional non-polarizable MD model;10 c Polarizable classical Drude oscillator model;10 d MDEC model, eqn (2.13), where emd and Eei are taken from b and c, respectively; e The relative error of Born factor, eqn (3.7). 2620 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011 View Online a. U 3.4 Dielectric properties Here we consider the dielectric properties of polar (alcohols) and non-polar (alkanes) solvents published recently.10'21 The dielectric constants of these solvents were calculated by using both the conventional non-polarizable model and polarizable classical Drude oscillator model.10'21 A significant improvement was obtained when the effects of electronic polarization were included, in particular for non-polar alkanes. The results of such calculations are reproduced in Tables 1 and 2. As far as the dielectric constant is concerned, MDEC model states that the static dielectric constant is simply a product of the value found in non-polarizable MD simulations and that of the electronic continuum, eqn (2.13). As non-polarizable simulations have already been done, we use the published results and check the above relation. The results are shown in Tables 1 and 2 for alcohols and alkanes, respectively, where MDEC model is compared with both polarizable Drude model and with experiment.50 To quantify the comparison, we introduce a parameter 0, which is a measure of how well a given dielectric constant reproduces the results of charging free energy calculations of spherical ions. Since the free energies are proportional to the corresponding Born factors, AG ~ (1 — 1/fi), the parameter 0 is defined as: 0: (l-l/fisim)-(l-l//xp) (1 - l/£exp) (3.7) In Tables 1 and 2, parameter 0 is shown for different types of simulations and for MDEC model. As predicted by the criterion in eqn (3.6) the traditional non-polarizable model satisfactorily reproduces the polarization effect in such polar media as alcohols (rmsd of 0 < 4%, Table 1); although, the polarizable classical Drude oscillator model21 and the MDEC model demonstrate much better agreement with the experiment (rmsd of 0 < 1% for both). In the case of non-polar media the traditional MD approach completely fails (rmsd of 0 ~ 100%, see Table 2); whereas, the polarizable Drude model10 and the MDEC model satisfactorily describe the polarization of neat alkanes (rmsd of 0 is 5.1 % and 4.3%, respectively). The MDEC approach appears to be even slightly favorable in this case. Thus, in the above examples MDEC approach performs quite well for both polar and non-polar media. We next explored the application of non-polarizable models to a low-dielectric interior of proteins, which have been studied in the past using standard MD simulations.12 For the dielectric permeability of the most internal region of cytochrome c, and several other proteins, Simonson et a/.12'15 reported a value around 1.5 (and even lower deeper inside). This value is apparently too low to be the actual dielectric constant of the protein; indeed it is lower than the pure electronic permeability £ei = 2 estimated for cytochrome c in the polarizable calculation.16 According to MDEC model, to obtain the total (static) dielectric constant, the results of non-polarizable simulations should be modified as given by eqn (2.13). Considering the value 1.5 as corresponding to emd, and using £ei = 2, for static dielectric constant we obtain the value £0 = 3.0, which is in agreement with the value 2.9 estimated by Muegge et al.51 In a related work, the non-polarizable simulations have been used for the analysis of dielectric properties of the interior of redox protein Cytochrome c oxidase52 (CcO). The charge insertion process has been studied that models deprotonation of His291 residue of CuB catalytic center in (dehydrated) CcO. The free energy data and corresponding dielectric properties are given in Table 3. The reaction-field energy obtained in the traditional MD technique was found to correspond to an un-physically low protein dielectric constant of 1.3. However, when the electronic (fiei = 2.0) polarization energy was added explicitly, as given by eqn (2.12), the microscopic reaction-field energy could be reproduced with a more realistic value of protein dielectric of 2.6. The estimated magnitude of the dry CcO dielectric constant in the region of active site is consistent with earlier results for cytochrome c which are 2.951 or 3.0 (the value15 corrected by eqn (2.13)). 3.5 Solvation free energy The comparison of the standard MD and MDEC simulations using GROMOS3 force field with experimental data for the electrostatic hydration free energy of polyatomic ions is given in Fig. 4. As seen, MDEC free energies reproduce experimental data within the experimental error which is typically about several kcal mol-1 for ions. As expected, the conventional non-polarizable MD also reasonably well reproduces the polarization effect in such a high dielectric media as liquid water. The quality of free energy simulations, however, is different in the low-dielectric media such as liquid cyclohexane or protein interior of protein Cytochrome c oxidase, which was described above, see Table 3. As predicted by the criterion in eqn (3.6), the traditional non-polarizable MD completely fails in the simulations of low-dielectric environment. The obtained in ref. 48 solvation free energies of Methyl and Propyl Table 3 Solvation free energy of ions in low-dielectric media obtained by different MD models System Methyl Guanidinium in cyclohexane48 Propyl Guanidinium in cyclohexane48 His291 in dehydrated Cytochrome c oxidase5 FF AGmc, kcal moF AG, kcal moP1 EMD Eel CHARMM -0.81 -27.5" 1.0C 2.0 2.050 CHARMM -0.75 -28.6" 1.0C 2.0 2.050 AMBER -15.7 -45.6fc 1.3C 2.0 2.6C ° Polarizable MD simulations. b MDEC, eqn (2.12). c Dielectric constant is estimated in continuum electrostatic calculations of the solvation free energy adjusting emd to reproduce AGMD or £0 to reproduce AG. This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2621 View Online a. U -60 u o CD < QM,D d/md JMD 9 B (RHFI QM„ ,/MDEC, k=0.9 Experiment IS s B t Me(CO/ MeNH3* Me2NH2* Me3NH+ NH," A MeO MeOH* J,_,_J_l boh; PhNH * CSH5NH* AG (exp), kcal/mol so/v Fig. 4 Charging free energy of polyatomic ions in aqua solution simulated in the ref. 28. Opened symbols correspond to the traditional MD simulations, while filled symbols stand for MDEC technique, eqn (2.10). The experimental values (dashed line) were obtained as a difference of total solvation energies for a given ion53 and its hydrocarbon counterpart.54 Partial charges and geometry of ions were found in self-consistent reaction-field quantum-chemical procedure.40 Circles correspond to AMI level of theory, while squares and triangles to RHF/6-31G**. Van der Waals radii of G43A force field3 scaled by the factor k — 0.9 (filled circles and squares) and k — 0.8 (filled triangles) were used to build the molecular cavity in MDEC computation of electronic free energy component. Guanidinium ions in cyclohexane are about zero (or of the order of thermal fluctuations ~kBT), i.e. the standard non-polarizable MD simulations result in no polarization effect at all, emd ~ 1, while free energies obtained48 in polarizable MD are several orders of magnitude higher reflecting the correct dielectric constant of cyclohexane e0 = 2.0.50 Thus, according to the criterion in eqn (3.6) and simulation data summarized above the standard non-polarizable MD technique significantly underestimate the solvation free energy of ions and medium dielectric constant in the low-dielectric materials. In contrast, MDEC simulations are correct for both high and low-dielectric media. 3.6 What are the microscopic interactions in the condensed phase? To test how well the macroscopic electronic continuum approximation of MDEC model describes a molecular solvent, with its microscopic structure and corresponding inhomogeneity of electronic density, we examined39 a model of two ions A~ and A + dissolved in benzene. In the non-polar solvent the nuclear component of polarization is negligible and the total effect is almost exclusively determined by the electronic polarization therefore the total potential of mean force (PMF) accurately represents the electronic screening effect. The low-dielectric environment is similar to that in the interior of a protein, or a lipid membrane. The solvent now is described by the polarizable Drude oscillator model,22 whereas ions are treated by a standard non-polarizable force field (Coulomb and Lennard-Jones interactions; the LJ parameters r,A 10 12 14 20 _ 0 O E -20 | -60 O. -80-I -100 Fig. 5 PMF for an ion pair A+ and A~ in benzene. The squares, circles and triangles stand for the results obtained with polarizable MD, non-polarizable CHARMM and CHARMM with scaled charges of the ions, respectively. Continuous curves are the least square fitting of the simulation points by the Coulomb function —l/er (with the Ewald correction, see ref. 39). For polarizable simulations (solid line), the effective dielectric constant £0 — 1.88; for non-polarizable simulations (dashed line), the dielectric constant emd — 1.16. The triangles correspond to non-polarizable CHARMM simulations with scaled charges by a factor 1/^(eq/emd) according to MDEC model. for ions correspond to those of Cl~ ion.) For such a system, we calculate the electrostatic part of the potential of mean force and compare the results with those of scaled and un-scaled CHARMM calculation, using the concepts of MDEC theory, see Fig. 5. The PMF gradient over r gives the average electrostatic force acting between charged particles A~ and A + in the bulk. The solvation free energy AG(r) of ions was evaluated by three alternative techniques: by polarizable MD, by the standard technique using non-polarizable CHARMM force field, and by using eqn (2.10) and CHARMM force field with scaled ion charges, according to MDEC model. As seen in Fig. 5, when the space between ionic spheres is larger than a size of solvent molecules, the effects of solvent microscopic structure becomes unimportant, and the average interaction, both in polarizable and non-polarizable models of benzene, can be approximated by a simple Coulomb law with an effective dielectric constant (obviously the LJ interactions are not important in this region). In case of polarizable Drude oscillator model for solvent benzene, the average interaction between ions is reproduced with an effective dielectric constant £0 = 1.88.** According to MDEC theory, eqn (2.13), the total dielectric constant of the medium e0 is a product of the electronic dielectric se\ (due to Drude-polarization of benzene molecules) and that of nuclei, emd. The latter was obtained in a separate simulation using non-polarizable CHARMM model of benzene; the corresponding value is emd = 1.16. According to eqn (2.13) then, the corresponding electronic dielectric constant of polarizable model of benzene is £el = e0/emd = 1.62. As seen in Fig. 5 (solid line), in perfect agreement with MDEC theory, the results of polarizable benzene simulations are reproduced by scaling charges of ions (by a factor 1 /vT~62) and running non-polarizable CHARMM simulations. Again, we see that all ** We notice that the experimental value of £0 for benzene is actually 2.3;50 the underestimated value of £0 is a consequence of the reduced polarizability parameter employed in the benzene model,22 which is ~20% lower than experimental benzene polarizability. This lack of parameterization, however, is not essential for the model test provided by the consistent choice of the underestimated value of eei. 2622 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011 View Online a. U effects of electronic polarization can be incorporated by scaling charges of ions with a factor 1 / The significant deviation of the results of standard non-polarizable MD from those of polarizable and MDEC techniques shown in Fig. 5, in fact, can be rationalized without the PMF simulations. Since scaling factor of each microscopic model is given by the corresponding dielectric constant (as defined in section 2.3), the PMF profiles are approximated by the corresponding Coulomb functions — \/e0r, —l/sMDr and — l/(eei£MD)r for the polarizable, non-polarizable CHARMM and MDEC techniques, respectively. Due to the relation (2.13), PMF functions for polarizable MD and MDEC should be the same, while, deviation from the CHARMM technique is estimated as (1 — l/fiei) • l/sMDr. Thus, for the low-dielectric media where £el = 2 and emd ~ 1, the deviation is ~ l/2r, which is significant even for larger separation distances (~16 kcal moP1 for r = 10 A). In the high-dielectric media (emd » 1), however, the difference will be much smaller. For instance, in water (fiei = 1.8, emd ~ 100 for TIP3P model55) the deviation will be just ~ 1 /225r, which is ~0.5 kcal moP1 even for the shortest separation r < 4A (contact ion pair: r < 2Rvdw). Since 0.5 kcal mol-1 is of the order of statistical uncertainty of MD the missing electronic screening effect is not noticeable in the standard non-polarizable simulations of water solutions. Thus, despite a complex nature of electronic polarization in a real system, the effect can be described reasonably well in different solvation conditions by a simple charge scaling procedure; this opens a way to modify the standard force fields so as to improve the description of their charged groups by effectively incorporating the electronic screening of charges. 3.7 Dynamics of salt bridges in proteins To demonstrate the significance of accounting for the electronic screening effect in protein dynamics simulations we modeled39 fluctuations of an important salt bridge (Arg438-PropD of heme a3, see the structure in ref. 39) in Cytochrome c Oxidase (CcO). This salt bridge (SB) controls water penetration to the hydrophobic (low-dielectric) cavity in the catalytic center of CcO.56'57 The strength of the electrostatic interaction of the salt bridge determines the rate of its opening/closing and, as a result, the probability of water transfer to/from the catalytic cavity. The distance d between 02D of A-propionate and 2HH2 of Arg43& has been chosen to characterize the fluctuations of the salt bridge gate during an MD run. The AMBER1 force field was used. The distribution functions for distance d obtained with scaled and original un-scaled charges are shown in Fig. 6a. Here no water in the cavity was included in the simulation. It is seen that the SB dynamics becomes qualitatively different once electrostatic interactions between the charged Arg43& + and the COO~ group of A-propionate are reduced by a factor of l/fiei (in the simulations, £ei = 2.0). In contrast to the standard MD simulations,57 the fluctuations observed in the scaled model are significantly larger, so that the internal water can now easily pass through the opened SB gate, and enter the catalytic cavity. In fact, during a 5 ns MD run with scaled charges, several such water transitions were observed. 41 3- ; \ a) 1- 3 . a 4 d, A b) d, A Fig. 6 Distribution functions of the distance d between 02D (A-propionate of heme a3) and 2HH2 (Arg438) of bovine CcO salt-bridge: (a) no water in the catalytic cavity; (b) 4 water molecules are added to the cavity. Dashed lines represent distributions obtained in the standard MD, while solid lines stand for the distributions obtained in the MD with scaled charges of the ionized groups. In Fig. 6b, the distribution functions of d are shown from simulations that included water in (and around) the catalytic cavity of the enzyme. As we already pointed out, the electronic screening affects not only charge-charge interactions, but interaction with water as well. Here TIP3P model is taken without modification; the charge scaling affects only the salt-bridge groups. As seen in Fig. 6b when the effects of electronic screening are included even more dramatic changes are observed. Thus, standard (un-scaled charges) MD simulations with and without water in the cavity lead to the conclusion that the salt bridge is formed 100% of the time; here stability of the salt bridge is quantified by the criterion d < 3 A, while, the bridge is observed only 98% or even 63% of the time in simulations with scaled charges without (see Fig. 6a) and with water (see Fig. 6b), respectively. It is clear, that the account for electronic screening of charged groups can give rise to qualitatively different results in simulations of proteins. As we have shown,39 this can be achieved in a computationally effective way by simple charge scaling of ionized groups in the protein. Unfortunately, there are no direct experimental data on the dynamics of the salt bridge discussed here to verify our proposal of electronic screening. However, as we argue in this paper such a scaling is obvious from theoretical point of view. An indirect comparison with an experiment, and support of charge scaling, is provided by some other computational studies, such as Zhu et a/.58 where a heuristic approximation for the charge scaling of ionized side chains (variable dielectric constant > 2) somewhat similar to ours was employed, which resulted in significant improvement in both side chain and loop prediction for protein conformations. This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2623 View Online a. U 4. Conclusions There is inconsistency in how the effects of electronic polariz-ability are treated in the commonly used non-polarizable empirical force fields.1"* The electronic screening effect, inherent for the condensed phase media, appears to be accounted for only for neutral moieties, whereas the charged residues are treated as if they were in vacuum. As a result, the electrostatic interactions between ionized groups are exaggerated by a factor of about 2. Also, an important electronic polarization term is typically neglected in simulations of solvation free energy or dielectric constant. The omission of the electronic contribution to the solvation energy is compensated by the exaggerated electrostatic interactions, but the complete compensation is possible only in high-dielectric media. The discussed here non-polarizable MDEC (Molecular Dynamics in Electronic Continuum) model provides a theoretical framework for systematic accounting of the effects of electronic polarization, and suggests a modification of the standard non-polarizable force fields1-4 to make them consistent with the idea of uniform electronic screening of partial atomic charges. In a few examples, we compared the traditional non-polarizable MD simulations with MDEC simulations, and demonstrated how the charges of ionized groups can be rescaled to correspond to MDEC model. The present theory states that the charges of ionized groups of the protein, as well as charges of ions, in simulations with existing non-polarizable potentials such as CHARMM, AMBER, etc should be scaled; i.e. reduced by a factor 1 /^/%{ (about 0.7), to reflect the electronic screening of the condensed medium relevant to biological applications. In the solvation free energy simulations the electronic part of the free energy (estimated by the continuum model) should be added explicitly to the nuclear part, which is obtained in non-polarizable MD simulations. The inclusion of electronic screening for charged moieties is shown to result in significant changes in protein dynamics and can give rise to new qualitative results compared with the traditional non-polarizable force fields simulations. Appendix Dielectric screening The reduction of electrostatic interactions in a dielectric medium with respect to that in vacuum is called the dielectric screening of interactions. The effective interaction between two localized groups of charges in a dielectric of e is defined via the potential of mean force (PMF): lfnt(r, e) = PMF = G(r, e) - G(r = oo, e), (Al) here r is an effective distance between two charge groups; G(r, e) is the total electrostatic free energy of the system; G(r = oo, e) is the sum of free energies of individual groups. The dielectric screening factor fir, e) then is defined as the ratio between interaction in dielectric and in vacuum: fr, e) = Uint(r, e)IÜn\r, e = 1). (A2) In general, a solution of the dielectric problem with appropriate boundary conditions is necessary to obtain the free energy G(r, e), effective interaction Vn\r, e) and screening factor fr, e). If molecules in the condensed medium are represented by the partial charges distributed inside the molecular cavity of general shape then the result obviously will depend not only on the dielectric property e and distance r but also on the molecular charge distributions p\(x), p2ix) and cavity shapes S\(5c), S2(x), where x stands for the Cartesian coordinates. The complexity of the problem is significantly reduced in the case of spherical molecular cavities. First, consider two point dipoles fli and jt2 at centers of the spherical cuts in the dielectric of e. The total free energy of the system is given by *¥(x; r, s)p(x)d x X- J ¥(x;r,£)P](x)d3x + l; I7! *¥(x; r, s)p2(x)d3x, (A3) where p(x) = p\(x) + p2ix) is the total external charge distribution and m(x; r, e) is the total potential induced by the charge distribution. In the limit of large distances, r » i?-radius of cavities, the higher-order re-polarization of the cavities in response to each other polarization can be neglected and the solution for the 2-sphere problem is expressed via solutions for single sphere problems: ^(x e V\, r, e) = (pin\x) + e V2; r, e) = R,e) = -e 3 \ 3 fa • r)fa ■ r) - fa ■ ^2)r2 2e+ 1 (A8) and the screening factor of electrostatic interactions defined by eqn (A2) is /(/>*,£)=£ 2e+ 1 (A9) The screening (A9), in fact, is true for arbitrary charge distributions p\{x), p2(x) in spherical cavities, for which the lowest multipole moment is the dipole. For a general distribution of charges qt the solution of the dielectric problem outside own cavity is given by the famous Kirkwood's expansion:60 N 00 2/+ 1 (Xi)' (/+ i) + /H+] iMcosö,) 2/+ 1 §£(/+!) + / 4i 1 ^ ,„ N — 2^ qimYim(6,(p) 2l+\r< (A10) here xt are the positions of charges qt in respect to the sphere center, 6t are the angles between r and xt vectors, Pt(z) and Y[m(6, tp) are the Legendre polynomials and spherical harmonics, respectively, and qlm = '}2qi(xi) Yfm(0i,q>ij are spherical multipole moments59 which are linearly related to the Cartesian multipole moments. In eqn (A10) the vacuum component of the potential is grouped by brackets and the prefactor reflects the screening of the field corresponding to each multipole. It is seen that at larger distances r » R the electric field is determined exclusively by moments qLm of the lowest nonvanishing multipole L and, therefore, the solution (A 10) for the arbitrary charge distribution is equivalent to the field from the point multipole qLm located at the center of the sphere. Thus, at larger distances r » R the screening factor between arbitrary charge distributions (with lowest moments L\ and L2) in two spheres is the same as that between point charges (if Lx = L2 = 0), point dipoles (if Lx = L2 = 1) or any appropriate lowest moments located at the center of these spheres. The screening factors between charged (L\ = 0) and dipolar (L2 = 1) as well as between charged (Lx = 0) and charged (L2 = 0) distributions are obtained in similar manner giving the following combination rules for the dielectric screening: fqq{r > R, e) = ^ (charge-charge fqd(r^R, e) fdd(r^> R,s)=s 2e- (charge—dipole) (All) 2e+ 1 (dipole—dipole) see though ref. 24 and 61. One should remember that the obtained scaling factors for charge-dipole or dipole-dipole interactions are only for specific models of molecular cavities—spheres, and for large distances, and therefore should not be taken directly as more appropriate than a straightforward factor 1/e, in particular because the shape of molecular cavity is ill-defined. However, it does show that the true scaling in some models can be different from a simple 1/e. Acknowledgements It is our pleasure to acknowledge many insightful discussions with Dr Marshall Newton. This work has been supported in part by the NSF grant PHY 0646273, and NIH GM054052. References 1 W. Cornell, P. Cieplak, C. Bayly, I. Gould, K. Merz, D. Ferguson, D. Spellmeyer, T. Fox, J. Caldwell and P. Kollman, /. Am. Chem. Soc, 1995, 117, 5179-5197; J. Wang, P. Cieplak and P. A. Kollman, /. Comput. Chem., 2000, 21, 1049-1074. 2 A. D. Mackereil, D. Bashford, M. Bellott, R. Dunbrack, J. Evanseck, M. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin and M. Karplus, /. Phys. Chem. B, 1998, 102, 3586-3616. 3 W. F. van Gunsteren, S. R. Billeter, A. A. Eising, P. H. Hunenberger, P. Kruger, A. E. Mark, W. R. P. Scott and I. G. Tironi, Biomolecular Simulation: The GROMOS96 Manual and User Guide, Vdf Hochschulverlag AG an der ETH Zurich, Zurich, Groningen, 1996. 4 W. L. Jorgensen and J. Tirado-Rives, /. Am. Chem. Soc, 1988, 110, 1657-1666; W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, /. Am. Chem. Soc, 1996, 118, 11225-11236. 5 H. Fröhlich, Theory of Dielectrics, Clarendon Press, Oxford, UK, 1949. 6 M. Karplus, Acc. Chem. Res., 2002, 35, 321-323. 7 G. Hummer, L. R. Pratt and A. E. Garcia, /. Phys. Chem., 1996, 100, 1206-1215. 8 M. V. Vener, I. V. Leontyev, Y. A. Dyakov, M. V. Basilevsky and M. D. Newton, /. Phys. Chem. B, 2002, 106, 13078-13088. 9 I. Vorobyov, V. M. Anisimov, S. Greene, R. M. Venable, A. Moser, R. W. Pastor and A. D. MacKerell, /. Chem. Theory Comput., 2007, 3, 1120-1133. 10 I. V. Vorobyov, V. M. Anisimov and A. D. MacKerell, /. Phys. Chem. B, 2005, 109, 18988-18999. 11 C. Oostenbrink, A. Villa, A. E. Mark and W. F. Van Gunsteren, /. Comput. Chem., 2004, 25, 1656-1676. 12 T. Simonson and C. L. Brooks, /. Am. Chem. Soc, 1996, 118, 8452-8458. 13 F. T. H. Leuwerink and W. J. Briels, /. Phys. Chem. B, 1997, 101, 1024-1034. This journal is © the Owner Societies 2011 Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 | 2625 View Online a. U 14 H. Nymeyer and H. X. Zhou, Biophys. J., 2008, 94, 1185-1193. 15 T. Simonson and D. Perahia, Proc. Natl. Acad. Sei. U. S. A., 1995, 92, 1082-1086. 16 T. Simonson and D. Perahia, /. Am. Chem. Soc, 1995,117,7987-8000. 17 T. A. Halgren and W. Damm, Curr. Opin. Struct. Biol., 2001, 11, 236-242. 18 S. W. Rick, S. J. Stuart and B. J. Berne, /. Chem. Phys., 1994, 101, 6141-6156. 19 S. J. Stuart and B. J. Berne, /. Phys. Chem., 1996, 100, 11934-11943; G. A. Kaminski, H. A. Stern, B. J. Berne and R. A. Friesner, /. Phys. Chem. A, 2004, 108, 621-627; S. Patel, A. D. Mackerell Jr. and C. L. Brooks III, /. Comput. Chem., 2004, 25, 1504-1514. 20 G. Lamoureux and B. Roux, /. Chem. Phys., 2003,119, 3025-3039. 21 V. M. Anisimov, I. V. Vorobyov, B. Roux and A. D. MacKerell Jr., /. Chem. Theory Comput., 2007, 3, 1927-1946. 22 P. E. M. Lopes, G. Lamoureux, B. Roux and A. D. MacKerell, /. Phys. Chem. B, 2007, 111, 2873-2885. 23 L. Krishtalik, A. Kuznetsov and E. Mertz, Proteins: Struct. Funct. and Bioinform., 1997, 28, 174-182. 24 L. I. Krishtalik, Chem. Phys., 2005, 319, 316-329. 25 T. W. Allen, O. S. Andersen and B. Roux, Biophys. Chem., 2006, 124, 251-267. 26 J. L. Gao and D. G. Truhlar, Annu. Rev. Phys. Chem., 2002, 53, 467-505; H. M. Senn and W. Thiel, Angew. Chem., Int. Ed., 2009, 48, 1198-1229. 27 I. V. Leontyev, M. V. Vener, I. V. Rostov, M. V. Basilevsky and M. D. Newton, /. Chem. Phys., 2003, 119, 8024-8037. 28 M. V. Vener, I. V. Leontyev and M. V. Basilevsky, /. Chem. Phys., 2003, 119, 8038-8046. 29 I. V. Leontyev and A. A. Stuchebrukhov, /. Chem. Phys., 2009, 130, 085102. 30 E. W. Castner and J. F. Wishart, /. Chem. Phys., 2010, 132; D. Bedrov, O. Borodin, Z. Li and G. D. Smith, /. Phys. Chem. B. 2010, 114, 4984^1997. 31 M. K. Gilson, A. Rashin, R. Fine and B. Honig, /. Mol. Biol., 1985, 184, 503-516; C. N. Schutz and A. Warshel, Proteins: Struct. Funct. and Bioinform., 2001, 44, 400^117. 32 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, /. Chem. Phys., 1983, 79, 926-935. 33 H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, /. Phys. Chem., 1987, 91, 6269-6271. 34 S. Miertus, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117-129. 35 D. P. Geerke and W. F. van Gunsteren, ChemPhysChem, 2006, 7, 671-678. 36 T. Simonson, Rep. Prog. Phys., 2003, 66, 737-787. 37 R. Platzman and J. Franck, Z. Phys., 1954, 138, 411^131. 38 M. Neumann and O. Steinhauser, Chem. Phys. Lett., 1984, 106, 563-569. 39 I. V. Leontyev and A. A. Stuchebrukhov, /. Chem. Theory Comput., 2010, 6, 1498-1508. 40 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. J. A. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Kiene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, Revision-B.O4, Gaussian, Inc., Pittsburgh PA, 2003. 41 L. Delle Site, R. M. Lynden-Bell and A. Alavi, /. Mol. Liq., 2002, 98-99, 79-86. 42 P. L. Silvestrelli and M. Parrinello, /. Chem. Phys., 1999, 111, 3572-3580; M. Sharma, R. Resta and R. Car, Phys. Rev. Lett., 2007, 98, 247401. 43 Y. S. Badyal, M. L. Saboungi, D. L. Price, S. D. Shastri, D. R. Haeffner and A. K. Soper, /. Chem. Phys., 2000, 112, 9206-9208. 44 J. G. Kirkwood, /. Chem. Phys., 1939, 7, 911-919; L. Onsager, /. Am. Chem. Soc, 1936, 58, 1486-1493. 45 I. V. Leontyev and A. A. Stuchebrukhov, /. Chem. Theory Comput., 2010, 6, 3153-3161. 46 F. J. Lovas, /. Phys. Chem. Ref. Data, 1978, 7, 1445-1750. 47 W. F. Murphy, /. Chem. Phys., 1977, 67, 5877-5882. 48 I. Vorobyov, L. Li and T. W. Allen, /. Phys. Chem. B, 2008, 112, 9588-9602. As given in the third column of Table 1, solvation free energies of ions in cyclohexane modeled by the standard MD are about several kcal mol-1; whereas, values obtained by the polarizable Drude model are about —30 kcal mol-1. Elimination of the non-electrostatic part from the results of non-polarizable MD gives the electrostatic solvation energies of —0.81 kcal mol-1 and —0.75 kcal moP1 for Methyl Guanidinium and Propyl Guanidinium, respectively. Later data were kindly provided by authors in a personal communication. 49 T. Simonson, J. Carlsson and D. Case, /. Am. Chem. Soc, 2004, 126, 4167^1180; G. Archontis and T. Simonson, Biophys. J., 2005, 88, 3888-3904. 50 CRC Handbook of Chemistry and Physics (Internet Version 2010), ed. D. R. Lide, CRC Press/Taylor and Francis, Boca Raton, FL, 2010. 51 I. Muegge, P. X. Qi, A. J. Wand, Z. T. Chu and A. Warshel, /. Phys. Chem. B, 1997, 101, 825-836. 52 I. V. Leontyev and A. A. Stuchebrukhov, /. Chem. Phys., 2009, 130, 085103. 53 J. Florian and A. Warshel, /. Phys.Chem. B, 1997, 101, 5583-5595. 54 J. M. Wang, W. Wang, S. H. Huo, M. Lee and P. A. Kollman, /. Phys. Chem. B, 2001, 105, 5055-5067. 55 P. Hochtl, S. Boresch, W. Bitomsky and O. Steinhauser, /. Chem. Phys., 1998, 109, 4927^1937. 56 D. M. Popovic and A. A. Stuchebrukhov, /. Am. Chem. Soc, 2004, 126, 1858-1871. 57 M. Wikstrom, C. Ribacka, M. Molin, L. Laakkonen, M. Verkhovsky and A. Puustinen, Proc. Natl. Acad. Sei. U. S. A., 2005, 102, 10478-10481. 58 K. Zhu, M. R. Shirts and R. A. Friesner, /. Chem. Theory Comput., 2007, 3, 2108-2119. 59 J. D. Jackson, Classical Electrodynamics, John Wiley & Sons, New York, 1998. 60 J. G. Kirkwood, /. Chem. Phys., 1934, 2, 351-361; S. Deng, /. Electrost., 2008, 66, 549-560. 61 V. M. Agranovich and M. D. Galanin, Electronic Excitation Energy Transfer in Condensed Matter, North-Holland, Amsterdam, 1982; R. S. Knox and H. van Amerongen, /. Phys. Chem. B, 2002, 106, 5289-5293; C. P. Hsu, G. R. Fleming, M. Head-Gordon and T. Head-Gordon, /. Chem. Phys., 2001, 114, 3065-3072; C. Curutchet, G. D. Scholes, B. Mennucci and R. Cammi, /. Phys. Chem. B, 2007, 111, 13253-13265. 2626 I Phys. Chem. Chem. Phys., 2011, 13, 2613-2626 This journal is © the Owner Societies 2011