Many-body effects and simulations of potassium channels

The electronic polarizability of an ion or a molecule is a measure of the relative tendency of its electron cloud to be distorted from its normal shape by an electric field. On the molecular scale, in a condensed phase, any species sits in an electric field due to its neighbours, and the resulting polarization is an important contribution to the total interaction energy. Electrostatic interactions are crucial for determining the majority of chemical–physical properties of the system and electronic polarization is a fundamental component of these interactions. Thus, polarization effects should be taken into account if accurate descriptions are desired. In classical computer simulations, the forces required to drive the system are typically based on interatomic interaction potentials derived in part from electronic structure calculations or from experimental data. Owing to the difficulties in including polarization effects in classical force fields, most of them are based just on pairwise additive interaction potentials. At present, major efforts are underway to develop polarizable interaction potentials for biomolecular simulations. In this review, various ways of introducing explicit polarizability into biomolecular models and force fields are reviewed, and the progress that might be achieved in applying such methods to study potassium channels is described.


Introduction
A force field refers to the functional form and parameter sets used to describe the potential energy of a system of particles.Force field functions and parameter sets are derived from both experimental work and high-level quantum-mechanical calculations.In particular, at the heart of the potential functions are the atomic partial charges and the Lennard-Jones radii, which surprisingly vary widely from one force field to another.The polarizability of a molecule or an ion is generally smaller in condensed phase than in gas phase because of the presence of a confining potential due to the surrounding species.The fact that the motion of every single atom influences and depends on the motion of the surrounding atoms, and therefore that coupled equations are required to describe the dynamics of the system, is referred as the many-body effect.Generally, most biomolecular force fields and simulations do not consider electronic polarization explicitly, using fixed permanent partial charges on the atoms or interaction sites.If polarization effects are included, it is usually in an empirical manner by a mean-field approach in which the gas-phase atomic charges are used to mimic the bulk-phase dipole moment in combination with a scaling factor reproducing many structural and dynamical properties.This represents a reasonable first-order approximation, leading to over-polarization when the ions or molecules are transferred to their condensed phases.However, in complex biological systems such as ion channels, where the ion chemical environments that should be described are very different, the mean-field approach might be inadequate.While a charge-charge interaction term coupled with Lennard-Jones forces would reproduce some of the effects present in such a system, they would not reproduce, for example, the changes in the electronic structure of carbonyls in the filter region of a potassium channel such as KcsA, which would be expected as a potassium ion proceeds through the selectivity filter region.This fact has been recognized in the development of potential functions specifically designed for ions occupying this area of the protein, which differ significantly from those for general usage in solutions (Roux & Berneche 2002).Owing to the approximations required in the creation of a classical force field, the optimal parametrization is the result of a compromise between, for example, an accurate representation of the microscopic energies and bulk solvation properties.For instance, despite its key role in biology, no universally satisfactory model potential has yet been built for liquid water (Yu & van Gunsteren 2005) even if several pairwise-additive, site-site models, such as TIP3P, TIP4P or SPC/E, have been described and extensively employed.

Potassium channels
Ion channels are a large and biomedically important family of integral membrane proteins (Hille 2001) that govern the electrical properties of the membranes of excitable cells such as neurons or muscle, although they are also found in the membranes of non-excitable cells and a wide range of organisms from viruses to plants.Ion channels form pores in biological membranes, and each time the channel opens, thousands of millions of ions diffuse down their electrochemical gradient across the membrane.K C channels are perhaps the most extensively studied family of ion channels, both experimentally and computationally.K C channels are selective to K C ions over other cations such as Na C .For instance, the voltage-gated KcsA channel has a K C flow of approximately 10 7 ions per second per channel versus only 10 4 Na C ions transported, even if the charge of these ions is the same, and they have very similar ionic radius (1.33A ˚for K C and 0.95 A ˚for Na C ).The selectivity filter of K C channels is the structural element essential to the permeation and the selectivity mechanisms of these membrane protein architectures (figure 1).Ions in the extracellular and intracellular milieux cross the cellular membrane through the protein and its selectivity filter.Four identical subunits that are symmetrically disposed forming the central pore are composed of a conserved signature sequence peptide, the TVGYG motif (Thr-Val-Gly-Tyr-Gly) (Morais-Cabral et al. 2001;Zhou et al. 2001).The carbonyl oxygens of the backbone of this motif point towards the lumen, orchestrating the movements of ions in and out of the channel.These key carbonyl oxygens together with the side-chain hydroxyl oxygen of the threonine residue define four equally spaced ion-binding sites.The selectivity filter is too narrow to fit an ion with its solvation shell.Thus, ions must be dehydrated to enter this region of the channel and thus K C ions replace the solvation shell by the carbonyl oxygens of the protein while travelling across the channel.Each of these protein sites binds the permeant cation with a tight-fitting cage of eight carbonyl oxygen atoms.This geometry resembles the solvation shell of hydrated K C ions.Selectivity of K C ions over other ions such as Na C was suggested to rely on the pore having a snug fit for only K C ions (Doyle et al. 1998;Gouaux & MacKinnon 2005).This suggestion would imply conformational rearrangements to discriminate two species of very close sizes.Later, it was proposed that K C permeation at the mouth level of the protein is facilitated owing to the fact that dehydration of K C ions is better compensated by reassociation with the protein carbonyls in contrast to Na C ions (Guidoni et al. 2000).More recently, it was suggested that the selectivity is determined by the ability of the channel to coordinate only some particular ion species and that the magnitude of the dipole moment of the carbonyl ligands is the predominant factor for the selectivity (Karreman & Eisenman 1962;Aqvist et al. 1992;Noskov et al. 2004).However, the dipole strength, of itself could lead to the smaller ions, Na C , being favoured in the pore over K C .Using molecular ) embedded in a lipid bilayer, with individual lipids rendered as cyan chains.The transmembrane a-helices are rendered as blue ribbons and the extracellular domain is coloured according to its structure (yellow, beta sheet; white, random coil; cyan, turns).K C ions are rendered as purple spheres, and water oxygen and hydrogen as red and white spheres, respectively.Only waters in the selectivity filter are shown.On the right, atomic detail of the selectivity filter (TVGYG motif ) in a representative conductive structure is presented.For visual clarity, atomic detail is given only of two diametrically opposed subunits.
dynamics simulations and ab initio calculations on small models, evidence was later presented suggesting that the magnitude of the dipole moment of the carbonyl ligands is not the predominant factor in K C selectivity, but that selectivity arises primarily from limitations on the way in which these dipoles can orient around an ion (Thomas et al. 2007;Varma & Rempe 2007;Varma et al. 2008).Based on quantum calculations, it was also proposed (Huetz et al. 2006) that a different loss of charge between K C and Na C across the channel could significantly contribute to this specific conduction and selectivity.Large polarization effects on the K C ions and the oxygen ligands were reflected in the values of the partial charges associated with these sets of atoms.These effects were larger for K C than for Na C , which meant that Na C ions are much more strongly trapped in the selectivity filter once they have managed to get in than K C ions.The weaker binding resulting for K C ions explained the higher translocation rate across the channel.

Approaches to polarization
Transferability of the functional form and parameters is an important feature of a force field, i.e. the same set of parameters should be capable of describing related systems.In the context of ion channel simulations, they usually do not allow for polarization effects (Roux & Berneche 2002), which is a possible limitation shared with a number of other studies of ion selectivity (Allen et al. 2000;Aqvist & Luzhkov 2000;Domene & Sansom 2003) and permeation (Berneche & Roux 2001;Domene et al. 2004;Holyoake et al. 2004;Furini et al. 2009).Certainly, electronic polarization will play a key role during the movement of the ions both in the aqueous solution and in the channel (Guidoni et al. 2000;Bucher et al. 2006) and it will merit further study.
A significant number of computational studies have been carried out on ion channels, potassium channels in particular, of which many have made use of classical force fields.These force fields, which include Chemistry at Harvard Macromolecular Mechanics (CHARMM; MacKerell et al. 1998), Assisted Model Building and Energy Refinement (AMBER; Cornell et al. 1995), GROningen MOlecular Simulation (GROMOS; Oostenbrink et al. 2004), Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA; Jorgensen & Tiradorives 1988) and many others, have a long track record of use in molecular modelling.However, the accuracy of such force fields for the simulation of ion channel systems has been called into question, for example by the work of Roux & Berneche (2002), who demonstrated the varying interaction energies between ions and N-methyl acetaldehyde with a range of force fields.While these authors are not entirely dismissive of the possibility for correctly calibrated molecular dynamics studies to produce meaningful results, the question is raised as to whether the use of models allowing for induced polarization might yield better data.Other authors have also raised the question of potential issues caused by the neglect of polarization (Tieleman et al. 2001;Domene & Sansom 2003;Jorgensen et al. 2007), and further studies demonstrating the importance of polarization in other ionic systems (Stuart & Berne 1996), as well as in ion channels themselves (Allen et al. 2004;Noskov et al. 2004;Bucher et al. 2006;Warshel et al. 2007), make this question all the more pertinent.As has been done in a number of cases, the use of ab initio methods avoids this difficulty (Guidoni et al. 2000;Ban et al. 2004;Bliznyuk & Rendell 2004;Compoint et al. 2004;Huetz et al. 2006;Kariev et al. 2007;Kariev & Green 2008).However, ab initio calculations on large molecular systems are still computationally expensive and somewhat impractical.In the past, a number of alternative methods have been suggested to incorporate polarization into such molecular models, a number of which attempt to model polarization by means of alterations to a standard classical framework.These methods that are illustrated in a qualitative manner in figure 2, are often divided into four broad approaches.
(a ) The fluctuating charge model Of these approaches, the fluctuating charge method of polarization, represented in figure 2b, was developed by Rick et al. (1994) on the basis of previous work by Rappe & Goddard (1991), is perhaps the closest to the original concept of a classical force field.Instead of applying the fixed charges typical of such a force field, charges are allowed to vary.Where standard classical force fields model the electrostatic energy of a system as the sum of atomic charge-charge interactions, the fluctuating charge method incorporates an additional electrostatic energy term, representing the energy required to modify the charge of an individual atom.The total electrostatic energy is then minimized with respect to the individual atomic charges, subject to a set of constraints, for example keeping constant the overall charges on specific molecules, or on the system as a whole.
Specifically, as expressed by Rappe and Goddard, in a system of N atoms with charges q 1 , q 2 , ., q N , the electrostatic energy is written as q a q b J ab ; ð3:1Þ where E a0 is the energy of atom a at zero charge; c 0 a is the electronegativity of atom a; and J ab is a term equal to the Coulomb interaction between unit charges at atom locations a and b if asb, and equal to the difference between the ionization potential and the electron affinity of the atom if aZb.Charges in this system are derived by calculating the N charges q 1 , q 2 , ., q N subject to the N constraints where Q is a fixed charge for the system as a whole.This model has found a significant number of applications to a variety of small-molecule systems, including models of water (Rick et al. 1994;Chen et al. 2000a,b;Rick 2001;Olano & Rick 2005), ethanol (Wang & Cann 2007) and ionic systems (Bryce et al. 1998;Ribeiro & Almeida 1999).More recently, the fluctuating charge model has been applied to classical force fields for modelling protein systems, being combined with the CHARMM force field (Patel et al. 2004), and in the work of Friesner and co-workers, being used in a stand-alone sense (Applequist et al. 1972;Banks et al. 1999), and in later work as part of a model including polarizable dipoles (Stern et al. 2001).
While the fluctuating charge model can be used as part of quantum mechanics/molecular mechanics (QM/MM) calculations (Field 1997) or in a stand-alone classical model (Olano & Rick 2005), the quantum mechanics/ molecular mechanics (induced-charge), QMMM(IC), model of polarization developed by Reynolds and co-workers (Ferenczy et al. 1997;Winn et al. 1997Winn et al. , 1999;;Ferenczy & Reynolds 2001;Illingworth et al. 2006) is explicitly designed for a QM/MM framework.Here, the charge distribution of the QM region is used to calculate induced dipoles on the MM points in the system, making use of the isotropic polarizabilities of Miller & Savchik (1979).If the field induced by the QM region is E a for some atom a, then the dipole vector m a is given by the equation where a a is the appropriate atomic polarizability.The effect of these dipoles is then approximated by modifying the surrounding point charges (Ferenczy 1991).This produces a new set of point charges, which are fed back into the QM/MM model, setting up an iteration procedure that converges in a few steps to give a polarized system.In common with the fluctuating charge model, polarization is here included without the addition of dipoles or additional off-atom charges, potentially allowing for easy incorporation into modelling software.This method has successfully been applied to energy calculations in both small-molecule (Illingworth et al. 2006) and enzymatic (Illingworth et al. 2008b) systems, and to the problem of ligand docking (Illingworth et al. 2008a).

(b ) The Drude oscillator model
The Drude oscillator model (Drude 1902), also referred to as the shell method (Dick & Overhauser 1958) or as the charge-on-spring method (Yu & van Gunsteren 2005), incorporates polarization into a classical force field by means of a set of additional massless charges attached, one to each of the atomic centres, by means of harmonic springs.Figure 2c gives a qualitative picture of this method.The charge at the atomic centre is correspondingly reduced by the value of the added charge, and the springs are kept short in order to recreate the effect of the original charge on the surrounding system.Spelling this out mathematically, suppose that a massless particle with a charge of KQ a is attached to the site of an atom a, by a short spring with force constant k a .The charge on the atom is correspondingly increased by Q a .When the spring is in equilibrium, the force due to the electric field will be equal to the force on the spring.Thus, where the displacement of the massless particle is denoted by the vector Dr a , E a Q a Z k a Dr a : ð3:5Þ The displacement of the massless particle creates a dipole, m a , in the direction of the displacement m a Z Q a Dr a : ð3:6Þ And combining equations (3.5) and (3.6), and rearranging term, gives Here the scalar term on the r.h.s. of the equation, equal to the square of the added charge divided by the spring constant, is analogous to the isotropic atomic polarizability of equation (3.4).
In the presence of an electric field, this addition of charged particles recreates the effect of induced dipoles at the atomic centres.As with the fluctuating charge model, this method has been applied in a significant body of work modelling water and other small-molecule systems (Sprik & Klein 1988;Saint-Martin et al. 2000;Hernandez-Cobos et al. 2005;Yu & van Gunsteren 2005;Yu et al. 2006;Anisimov et al. 2007).Considering the role of polarization in ion channels, Roux and co-workers have developed a model for water on this basis (Lamoureux et al. 2003(Lamoureux et al. , 2006)), and applied it to a model of the hydration of K C ions (Whitfield et al. 2007).The Drude oscillator model was found to give results consistent with ab initio calculations, and in good agreement with the experimental measurements of the radial distribution function, thereby suggesting the potential of this method for application to systems of greater complexity.In addition to studies of small-molecule systems, polarization models based on the Drude oscillator model have also been applied to protein systems, being implemented in a range of common classical molecular dynamics frameworks (Lindahl et al. 2001;Lamoureux & Roux 2003;Yu & van Gunsteren 2004;Anisimov et al. 2005).

(c ) Induced dipoles
A third model for polarization in classical systems incorporates explicit dipoles at each of the atomic points, illustrated by figure 2d.The dipole m a at a point a with polarizability a a is expressed as where T ab is a tensor calculated from the relative positions of the atoms a and b.Note again, the similarity between this equation and equation (3.4), and the additional term modelling the dipole-dipole interaction.The polarizability, a a , is not necessarily isotropic, and therefore takes the form of a matrix.
Equation (3.8) is solved by rearranging into the form where the components T ab and a K1 a are 3!3 matrices and the m a and E a terms are 3!1 vectors.This gives a system of 3N equations, solved for the dipole terms m a , by inverting the 3N!3N matrix.
Early work employing this method was carried out by Applequist et al. (1972) and Warshel & Levitt (1976).Thole (1981), incorporated a damping scheme to the dipole-dipole interactions in order to avoid a situation in which dipoles placed too close together lead to an infinite polarizability.This work in turn was the basis for a series of publications including those by Ponder and co-workers (Ren & Ponder 2002, 2003;Schnieders et al. 2007), in which induced dipoles were used to study small-molecule systems and solvation effects.Indeed, the method of modelling polarization by explicit induced dipoles has been applied to models of a wide range of systems (Borodin & Smith 2006a,b;Mayer et al. 2006), including force fields designed for protein systems (Wang et al. 2006).Calculations by Baucom et al. (2004) have demonstrated that this brings about a meaningful improvement in the modelling of a DNA strand.
A more recent adaptation of this method replaces the framework of point charges and dipoles with Gaussian functions (Chelli & Procacci 2002;Paricaud et al. 2005), for example, through the use of s-orbital functions in place of point charges, and p-orbital functions in place of the point dipoles.This offers an alternative to Thole's damping method, which avoids the polarization catastrophe so long as the Gaussian functions are sufficiently diffuse.Use of such Gaussian functions has been shown to give an improved fit to molecular polarizability tensors when compared with Thole's model (Elking et al. 2007).
Also along these lines, the sum of interactions between fragments computed ab initio (SIBFA) method, developed by Gresh et al. (1984), is not strictly a classical method, in that it incorporated more than just polarization into the classical model, including, for example, an electrostatic term based upon ab initio calculations on fragments of the molecules of a system.However, it can be broadly fitted into this category, in that it makes use of point dipoles to model polarization.The SIBFA method has been used in a number of studies of ligand docking (Antony et al. 2005;Miller Jenkins et al. 2007), as well as in other applications, summarized in a recent review paper (Gresh et al. 2007).
The incorporation of atom-centred dipoles into a classical force field is an approximation to the higher order representation of the electrostatics of a molecule that can be achieved by means of multipole series.At the most basic level, the electrostatics of a small molecule can be represented by a single point with fixed charge, for example q Z X N aZ1 q a ; ð3:11Þ where q a is the charge on atom a, in the molecule.
A more detailed description is given by adding to this a dipole term aZ1 q a a i ; i Z 1; 2; 3; ð3:12Þ where the atom a has coordinates (a 1 , a 2 , a 3 ).
Higher order terms can be added to the description, the next being the quadrupole term where d ij is the Kronecker delta function.Higher order terms can also be calculated, the sum to infinity of all the terms giving a complete description of the electrostatics.While potentially useful, this approach suffers from slow convergence, a large number of terms being required to give a good description of the behaviour of a molecule.A significant development, therefore, in the use of multipoles, was the work of Stone on the derivation of atom-centred multipoles from a charge distribution (Stone 1981).Using a set of multipole series, centred on the atoms of a molecule, led to much faster convergence, and thereby a more practical tool for calculation.This was later extended to model the response of a molecule to an external electric field (Stone 1985).More recent development of this work has included the development of distributed multipole methods for evaluating the induction energies in condensed phase of small organic molecules, a project with direct application to organic crystal structure prediction (Misquitta & Stone 2008a,b;Misquitta et al. 2008;Welch et al. 2008).

(d ) QM-based polarization
A fourth method for including polarization, distinct from those above, involves the use of ab initio calculations on small parts of the molecules concerned, aiming to retain the accuracy of the full QM model, without carrying out an ab initio calculation on the whole system.Whereas the methods above broadly take the approach of adding components to a classical force field, this method takes a different ideological direction, beginning with a QM model that is simplified to be computationally tractable.An early example of this concept is the work of Field (1997) who incorporates semi-empirical molecular orbital methods into a QM/MM framework, applying the method of flexible charges to the classically modelled atoms.This approach has been taken in the MODEL potential method of Gao and co-workers (Gao 1997), in which individual molecules of a system are described at a semi-empirical level, with interactions between molecules being calculated on a QM/MM basis.This was used to develop a model for water (Gao 1998), and more recently to describe a force field suitable for calculations on peptides (Xie & Gao 2007).
In another approach, a fit is performed on a representative training set of molecules to reproduce some of their electronic properties calculated by accurate ab initio methods.Then the resulting atomic parameters are transferable both to new isolated molecular systems and to complexes in condensed-phase simulations (Kaminski et al. 2004;Chelli et al. 2005).A more accurate quantum-mechanical modelling of electronic polarization would be prohibitively expensive.
Simulations of biological systems often necessitate the use of large simulation cells and run times, which puts them outside the range of direct ab initio simulations.Potentials may then be parametrized by fitting the predicted forces and multipoles to a large body of information generated from ab initio calculations.The resulting potentials are predictive of ab initio accuracy and have a high degree of transferability (Madden et al. 2006).Compared with the almost parameter-free first-principles molecular dynamics in which electrons are explicitly represented, larger simulation cells and longer simulation times can be reached.A drawback although is that additional parameters are introduced in the interaction potential, and it is not always straightforward to choose them in a physically meaningful way.In the case of the polarizability, gas-phase values have often been measured experimentally, but values in liquid phase remain undetermined.

Discussion
Despite a considerable body of work on the problem, consensus has not yet been reached on a universally applicable model for polarization (Jorgensen 2007).Each of the classical approaches to modelling polarization is open to criticism.For example, despite recent efforts (Soteras et al. 2007), there is, as yet, no one agreed-upon method of defining atomic polarizabilities, essential for a dipolebased model.Drude oscillator models likewise present a question as to how to generate a correct parametrization.While if point charges alone are used to recreate polarization effects, there exists a geometrical restriction in polarization inherent to the location of the points.Although the use of additional off-atom charges provides a potential solution to this, this itself presents an issue of parametrization (Illingworth et al. 2006).QM-based models perhaps provide a better alternative, although they do not quite fit the ideal of incorporating polarization into a classical system.
Given the diversity of approaches summarized above, it is difficult to suggest which among them might prove the most useful for specific application to ion channels.In classical simulations where polarization is not included, there exist a number of independently derived force fields that have each produced valuable results, and, similarly, different models of polarization might each give insight into these systems.With regard to ion channel calculations, a simple modification of the ionic charge, or the addition of a charge-on-spring or dipole term to the ion itself may not be complex enough to capture the subtlety of the protein-ion and protein-water interactions, and therefore, for the ions at least, a QM-based model may be required to get the best result.If so, this would perhaps favour the use of QM-based or fluctuating charge-type approaches, these being more straightforwardly compatible with ab initio or QM/MM calculations.However, understanding of how best to incorporate polarization into models of these systems will most probably be gained through future experimentation.
Here we note that, while no one method of incorporating polarization into classical models has been agreed upon, there is significant agreement on the importance of polarization in ion channels, as noted above.While polarizable classical force fields have been applied to a great number of small-molecule systems, and to a variety of protein simulations, the number of applications to ion channels is relatively small.Aqvist & Warshel (1989) implemented a method of Langevin dipoles in order to model the solvation of sodium ions by gramicidin A. Later studies have adopted similar methods in order to model polarizability in the KcsA channel (Luzhkov & Aqvist 2000;Burykin et al. 2003).Bucher et al. (2006) showed that there is a significant transfer of charge from the backbone atoms of the KcsA potassium channel, to the bound cations, thus the lack of polarization is a clear source of error.Given the recognized importance of polarization in these systems, ion channels perhaps provide a fertile ground for the future development and application of classical polarization techniques.
Various ways of introducing explicit polarizability into biomolecular models and force fields have been reviewed and discussed.They undoubtedly represent a significant step forward in devising new and physically realistic force fields for simulations of complex systems.However, it remains the case that all available polarizable force fields are based on classical electrostatics, none of which explicitly consider important non-classical effects such as polarization-exchange coupling or charge transfer.
Transferability of the functional form and parameters is an important feature of a force field; that is, the same set of parameters should be capable of describing related systems.In the context of ion channel simulations, they usually do not allow for polarization effects (Roux & Berneche 2002), which is a possible limitation shared with a number of other studies of ion selectivity (Allen et al. 2000;Aqvist & Luzhkov 2000) and permeation (Berneche & Roux 2001).In the absence of a fully developed polarizable force field for ions, water, protein and lipids, the studies carried up to now have delineated some of the main features of ion channels.However, it is expected that future research directed at using polarizable force fields and QM/MM and ab initio simulation methods, with the advent of linear-scaling quantum-mechanical methods, will lead to further insights into these biologically important systems.
C.D. thanks The Royal Society for a University Research Fellowship.This work was supported by grants from The Leverhulme Trust and the EPSRC.

Figure 2 .
Figure 2. Qualitative description of approaches to modelling polarization.(a) Neglect.Classical modelling of systems is regularly carried out with fixed atomic charges.(b) Fluctuating charge model.Atomic charges are modified to reproduce the effect of polarization.(c) Drude oscillator model.A point charge is attached to each atom by means of a spring, and moved, in order to reproduce polarization effects.(d ) Induced dipole model.An explicit dipole term is added to each point charge.(e) Quantum-mechanical model.A wave function is calculated in order to model polarization on an electronic level.